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Abstract: In this article, we analyze the stability of various numerical schemes for differential 
models of viscoelastic fluids. More precisely, we consider the prototypical Oldroyd-B model, for 
which a free energy dissipation holds, and we show under which assumptions such a dissipation 
is also satisfied for the numerical scheme. Among the numerical schemes we analyze, we consider 
some discretizations based on the log-formulation of the Oldroyd-B system proposed by Fattal and 
Kupferman in |15!, which have been reported to be numerically more stable than discretizations 
of the usual formulation in some benchmark problems [23j . Our analysis gives some tracks to 
understand these numerical observations. 
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Schemas dissipatifs en energie libre pour le modele 

d'Oldroyd-B 

Resume : Dans cet article, nous analysons la stabilite de divers shemas numeriques pour des 
modeles differentiels de fluides viscoelastiques. Plus precisement, on considere le modele d'Oldroyd- 
B comme prototype pour lequel on sait montrer qu'une energie libre est dissipee, et on verifie 
sous quelles hypotheses le shema numerique satisfait une propriete semblable. Parmi les shemas 
numeriques analyses figurent des discretisations fondees sur la reformulation logarithmique du 
systeme Oldroyd-B telle que proposee par Fattal et Kupferman dans [T5], reformulation qui a 
permis d'obtenir des resultats numeriques apparemment plus stables que les formulations usuelles 
dans certains cas tests [23] ■ Notre analyse donnc ainsi des pistes pour la comprehension de ces 
observations numeriques. 

Mots-cles : Fluides viscoelastiques ; Nombre de Weissenberg ; Stabilite ; Entropie ; Methodes 
des elements finis ; Methode de Galerkin discontinue ; Methode des caracteristiques 
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1 Introduction 

1.1 The stability issue in numerical schemes for viscoelastic fluids 

An abundant literature has been discussing problems of stability for numerical schemes discretizing 
models for viscoelastic fluids for over twenty years (see [27l EHl [Ml [32] for a small sample). Indeed, 
the numerical schemes for macroscopic constitutive equations are known to suffer from instabilities 
when a parameter, the Weissenberg number, increases. Such an unstable behavior for viscoelastic 
fluids at moderately high Weissenberg (Wi < 10) is known to be unphysical. 

Many possible reasons of that so-called high Weissenberg number problem (HWNP) have been 
identified [IHl [33 [3S1 [S] • However, these results have not led yet to a complete understanding of 
the numerical instabilities, despite some progress [Tni[53]. Roughly speaking, we can distinguish 
between three possible causes of the HWNP: 

1. Absence of stationary state: In many situations (flow past a cylinder, 4:1 contraction), the 
existence of a stationary state for viscoelastic models is still under investigation. It may 
happen that the non-convergence of the numerical scheme is simply due to the fact that, for 
the model under consideration, there exists no stationary state. 

2. Bad model: More generally, the instabilities observed for the numerical scheme may originate 
at the continuous level, if the solution to the problem indeed blows up in finite time, or if it 
is not sufficiently regular to be well approximated in the discretization spaces. 

3. Bad numerical scheme: It may also happen that the problem at the continuous level indeed 
admits a regular solution, and the instabilities are only due to the discretization method. 

In this paper, we focus on the third origin of instabilities, and we propose a criterion to test the 
stability of numerical schemes. More precisely, we look under which conditions a numerical scheme 
does not bring spurious free energy in the system. We concentrate on the Oldroyd-B model, for 
which a free energy dissipation is known to hold at the continuous level (see Theorem 12.11 below 
and [2^) and we try to obtain a similar dissipation at the discrete level. It is indeed particularly 
important that no spurious free energy be brought to the system in the long-time computations, 
which are often used as a way to obtain the stationary state. 

The Oldroyd-B system of equations is definitely not a good physical model for dilute polymer 
fluids. In particular, it can be derived from a kinetic theory, with dumbbells modeling polymer 
molecules that are unphysically assumed to be infinitely extensible. But from the mathematical 
viewpoint, it is nevertheless a good first step into the study of macroscopic constitutive equations 
for viscoelastic fluids. Indeed, it already contains mathematical difficulties common to most of the 
viscoelastic models, while its strict equivalence with a kinetic model allows for a deep understanding 
of this set of equations. Let us also emphasize that the free energy dissipation we use and the 
numerical schemes we consider are not restricted to the Oldroyd-B model: they can be generalized 
to many other models (like FENE-P for instance, see [10]), so that we believe that our analysis 
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can be used as a guideline to derive "good" numerical schemes for many macroscopic models for 
viscoelastic fluids. 



1.2 Mathematical setting of the problem 

We consider the Oldroyd-B model for dilute polymeric fluids in d-dimensional flows {d = 2,3). 
Confined to an open bounded domain T) C M"^, the fluid is governed by the following nondimcn- 
sionalized system of equations: 

f du \ 
Re ( — + M • Vu ] = -Vp + (1 - e) Am + div r, 



divM = 0, (1) 



d 1 

^ + (M . V)r = (Vm)x + T{Vuf - + 4 



Vm + VuF 



where u : {t,x) £ [0,T) x X> ^ u{t,x) e M'' is the velocity of the fluid, p : {t,x) e [0,T) x V 
p{t, a;) e M is the pressure and t : {t, x) e [0,T) x T) T{t, x) E M'^^'* is the extra-stress tensor. 
The following parameters are dimensionless: the Reynolds number Rc G M+, the Weissenberg 
number Wi G M*^ and the elastic viscosity to total viscosity fraction e G (0, 1). 

In all what follows, we assume for the sake of simplicity that the system (jlj is supplied with 
homogeneous Dirichlct boundary conditions for the velocity m: 

M = on dV. (2) 

Therefore, we study the energy dissipation of the equations ([T]) as time goes, that is, the way (m, t) 
converges to the stationary state (0,0) (equilibrium) in the longtime limit t oo. For free energy 
estimates under non-zero boundary conditions, we refer to |24j . 

Local-in-time existence results for the above problem have been proved in the bounded domain 
[0,T) X V when the system is supplied with sufficiently smooth initial conditions u{t = 0) and 
x(t — 0) (see [T^ and [T7j for instance). Moreover, global-in-time smooth solutions of the system ((ij 
are known to converge exponentially fast to equilibrium in the sense defined in [24j . Let us also 
mention the work of F.-H. Lin, C. Liu and P.W. Zhang [34] where, for Oldroyd-like models, local-in- 
time existence and uniqueness results are proven, but also global-in-time existence and uniqueness 
results for small data. Notice that more general global-in-time results have been collected only 
for a mollified version of the Oldroyd-B system ([T]) (see [1]), for another system close to ((TJ, the 
co-rotational Oldroyd-B system (see [35]), or in the form of a Beale-Kato-Majda criterion when 
P = (see US])- Even though the question of the global-in-time existence for some solutions of 
the Oldroyd-B system ([T]) is still out-of-reach, it is possible to analyze global-in-time existence for 
solutions to discretizations of that system. This will be one of the output of this article. 
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1.3 Outline of the paper and results 

We will show that it is possible to build numerical schemes discretizing the Oldroyd-B system ([T])- 
such that solutions to those discretizations satisfy a free energy estimate similar to that estab- 
lished in \24\ I20j for smooth solutions to the continuous equations. Our approach bears similarity 
with [36] , where the authors also derive a discretization that preserves an energy estimate satisfied 
at the continuous level, and with [32j . where another discretization is proposed for the same energy 
estimate as in [31]. Yet, unlike the estimates in El], our estimate, the so-called free energy 
estimate derived in [211 |2Qj, ensures the long-time stability of solutions. As mentioned above, 
long-time computations are indeed often used to obtain a stationary state. 

We also analyze discretizations of the log- formulation presented in [T51 (TB] , where the authors 
suggest to rewrite the set of equations ((ij after mapping the (symmetric positive definite) confor- 
mation tensor: 

Wi 

a^I+—T (3) 

e 

to its matrix logarithm: 

■0 = In CT. 

In the following, we assume that: 

cr(t — 0) is symmetric positive definite, (4) 

and it can be shown that this property is propagated in time (see Lemma 12.11 below) , so that 
tp is indeed well defined. The log-formulation insures, by construction, that the conformation 
tensor always remains symmetric positive definite, even after discretization. This is not only 
an important physical characteristic of the Oldroyd-B model but also an essential feature in the 
free energy estimates derived beneath. Besides, this log-formulation has indeed been observed to 
improve the stability of some computations [SO] [161 ES] . It is thus interesting to investigate whether 
the numerical success of this log-formulation may be related to the free energy dissipation. 
The main outputs of this work are: 

• One crucial feature of the numerical scheme to obtain free energy estimates is the discretiza- 
tion of the advection term {u ■ V)t (or {u ■ V)i/') in the equation on the extra-stress tensor. 
We will analyze below two types of discretization: characteristic method, and discontinuous 
Galerkin method (see Sections [H and O . 

• To obtain free energy estimates, we will need the extra-stress tensor to be discretized in a 
(elementwise) discontinuous finite element space (see Sections |4] and [5]). Yet, using a colloca- 
tion method to compute bilinear forms, and with special care dedicated to the discretization 
of the advection terms, it would still be possible to use continuous finite element spaces, as 
it will be shown in a future work [3]. 

• The existence of a solution to the numerical schemes that satisfy a free energy estimate will 
be proved whatever the time step for the log- formulation in terms of if), while it will be shown 
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under a CFL-like condition for the usual formulation in terms of t (sec Section ^ . This is 
due to the necessary positivity of the conformation tensor cr for the free energy estimates in 
the usual formulation to hold, and may be related to the fact that the log-formulation has 
been reported to be more stable than the formulation in terms of r (see [23]), though the 
uniqueness of solutions is still not ensured in the long time limit and bifurcations may still 
occur. This result will also be re-investigated in the future work [5] for a slightly modified 
formulation of the Oldroyd-B system of equations in terms of r, somewhat closer to the usual 
formulation than the log-formulation. 

Notice that we concentrate on stability issues. All the schemes we analyze are of course consistent, 
but we do not study the order of consistency of these schemes, neither the order of convergence. 
Again, this will be included in the future work [3] for the slightly modified formulation of the 
Oldroyd-B system of equations in terms of r mentioned above. 

Let us now make precise how the paper is organized. In Section [21 we formally derive the free 
energy estimates for the Oldroyd-B set of equations and for its logarithm formulation, in the spirit 
of [20] . Then, Section [3| is devoted to the presentation of a finite element scheme (using piecewise 
constant approximations of the conformation tensor and its log-formulation, and Scott- Vogelius 
finite elements for the velocity and pressure) , that is shown to satisfy a discrete free energy estimate 
in Section m Some variants of this discretization are also studied, still for piecewise constant stress 
tensor, and a summary of the requirements on the discretizations to satisfy a free energy estimate 
is provided in Tables [1] and [H We show in Section [5| how to use an interpolation operator so as to 
adapt the previous results with piecewise linear approximations of the conformation tensor and its 
log-formulation. Finally, in Section [H we show how the previous stability results can be used to 
prove long-time existence results for the discrete solutions. Some numerical studies are needed to 
illustrate this numerical analysis, and this is a work in progress. 

1.4 Notation and auxiliary results 

In the following, we will make use of the usual notation: L'^i'D) — {f : V W, J |/p < oo}, 



H\V) = {f:V-.Rj + |V/|2 < ex.}, H^V) = {/ : P ^ Rj\f\' + \Vf\^ + IWP < oo}. 



C([0,T)) for continuous functions on [0,r) and C^{[Q,T)) for continuously differentiable functions 



We will denote by r : cr the double contraction between rank- two tensors (matrices) x, cr G 



Notice that if r is antisymmetric and cr symmetric, r : cr = 0. 

The logarithm of a positive definite diagonal matrix is a diagonal matrix with, on its diagonal, 
the logarithm of each entry. We define the logarithm of any symmetric positive definite matrix cr 
using a diagonal decomposition cr ~ R^AR of cr with R an orthogonal matrix and A a positive 



on [0,T). 




l<i,j<d 
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definite diagonal matrix: 

\n(T = R^\nAR. (5) 

Although the diagonal decomposition of cr is not unique, ([5]) uniquely defines In cr. The matrix log- 
arithm bijectively maps the set of symmetric positive definite matrices with real entries iS^(M''^'*) 
to the vector subspace S{M.'^^'^) of symmetric real matrices, where it is exactly the reciprocal 
function of the matrix exponential. 

We will make use of the following simple algebraic formulae, which are proved in Appendix lA.il 
and[0] 

Lemma 1.1. Let cr and t be two symmetric positive definite matrices. We have: 

tr In cr = In dot cr, (6) 

cr — In cr — / is symmetric positive semidefinite and thus tr(cr — In cr — 7) > 0, (7) 

cr + — 21 is symmetric positive semidefinite and thus tr(cr + cr^ — 21) > 0, (8) 

tr(crr) = tr(Tcr) > 0, (9) 

tr ((cr - t)t-^) = tr(crx"i - I) > lndet(crx"i) = tr (In cr - Inr) , (10) 

tr ((Incr - Inr) cr) >tr(cr- r). (11) 
We will also use the usual Jacobi's formula: 
Lemma 1.2. For any symmetric positive definite matrix cr(t) £ ([0,T)), we have \ft G [0, T); 

— Incr : cr = tr cr— Incr = — tr cr. (13) 
dt J \ dt ) dt ^ ' 

2 Formal free energy estimates at the continuous level 

We are going to derive free energy estimates for two formulations of the Oldroyd-B system in 
Theorems 12.11 and 12.21 An important corollary to these theorems is the exponential convergence 
of the solutions to equilibrium in the longtime limit. In all that follows, we assume that (u,p, r) 
is a sufficiently smooth solution of problem ([T]) so that all the subsequent computations are valid. 
For example, the following regularity is sufficient: 

(«,p,T)e (Ci([0,T),i/2(p)))'^x (C0([0,r),i7i(2?))) X (Ci([0,T), 6-1(1?)))'"', (14) 
where we denote, for instance by ((7^(1?))' a vector field of dimension d with C^{T>) components. 
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2.1 Free energy estimate for the Oldroyd-B system 
2.1.1 Conformation-tensor formulation of the Oldroyd-B system 

Recall that the conformation tensor cr is defined from the extra-stress tensor t through the following 
bijective mapping: 

Wi ^ ' 

With this mapping, it is straightforward to bijectively map the solutions of system (jlj with those 
of the following system: 

Re — + M • Vtt = -Vp+ (1 - e)AM+ — div<T, 
\ ot J Wi 

divM = 0, (15) 

d(T 1 
— + (m • V)cr = {Vu)n- + (T{Vuf - - I)- 



Notice that with such an affine mapping, the solution cr to system (|15|) has the same regularity 
than T solution to system ^ , which is that assumed in ([H]) for the following manipulations. 

2.1.2 A free energy estimate 

Let us first recall a free energy estimate derived in [2H [2U] . The free energy of the fiuid is defined 
as the sum of two terms as follows: 

n«,^) = f /j«P + ^/^tr(.-ln.-/). (16) 

The kinetic term / jup is always non negative. Besides, we have the following lemma (see Ap- 
J-D 

pendix iBl or [22] for a proof): 

Lemma 2.1. Let cr G (C^ ([0, T), C^(I?)))''^'' he a smooth solution to the system (fT5| . Then, if 
the initial condition cr(t = 0) is symmetric positive definite (everywhere in V), the solution cr(t) 
remains so at all times t G [0,T) and for all x CzV. In particular, the matrix cr{t) is invertible. 

From Lemma |2. II and the equation ([7]), the entropic term / tr((T — In cr — /) is thus well defined 

Jv 

and non negative, provided cr{t — 0) is symmetric positive definite. 

The free energy is an interesting quantity to characterize the long-time asymptotics of the 
solutions, and thus the stability of system (jlSp . A priori estimates using the free energy are 
presented in Pl] for micro-macro models (such as the Hookean or the FENE dumbbell models) 
and in [20] for macroscopic models (such as the Oldroyd-B or the FENE-P models). Similar 
considerations can be found in the physics literature about thermodynamic theory for viscoelastic 
models (see (31 113 [H] ) . 

For the sake of consistency, we recall results from [20] : 
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Theorem 2.1. Let {u,p, cr) he a smooth solution to system (|15p supplied with homogeneous Dirich- 
let boundary conditions for u, and with symmetric positive definite initial condition (T(t = 0). The 
free energy satisfies: 



(17) 



IV 

From this estimate, we get that F{u, cr) decreases exponentially fast in time to zero. 

Proof of Theorem (|2.ip . Let (UjP, c) be a smooth solution to system (jlSp . with symmetric positive 
definite initial condition a(t = 0). We first compute the inner product of the Navier-Stokes equation 
with the velocity: 

Then, taking the trace of the evolution equation for the conformation tensor, we obtain: 

d 



dt 



/ trcr = 2 / Vu :<T ^ / tr((T-/). 

Jv Jv Wi7i5 



(19) 



Last, remember that smooth solutions cr are invertible matrices (Lemma I2.ip . Thus, contracting 
the evolution equation for cr with cr^^, we get: 



Using dEl) with (T eC^ {Vx [0,r),5^(M''^'^)), we find: 

which can be combined with ([^0]) to get, using tr(VM) = divit = and m = on dT>: 



(20) 



We now combine ([18]) +2Wi^ CS) -2Wi^ (EJ) to obtain ((T7|) 



(T-i-2/) = 0. 



Since, by (O, we have tr(cr + cr^^ — 27) > 0, then F{u, cr) decreases in time. Moreover, by ([7|) 
applied to cr~^, we have tr(cr — In cr — /) < tr(cr + cr~^ — 27). So, using the Poincare inequality 
which states that there exists a constant Cp depending only on T> such that, for all u e Hq{T>), 



V Jv 



2 



we finally obtain that F{u, cr) goes exponentially fast to 0. Indeed, we have from (fT7|l : 
iLf(„,.)<-i_£|^|„p__^|^ + 2/), 
. /2(l-e) 1 \ 
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so that, by a direct application of Gronwall's lemma, we get: 

F{u, a) < F{u{t = 0), a{t = 0)) exp (mm I * 



□ 



2.2 Free energy estimate for the log-formulation of the Oldroyd-B sys- 
tem 

2.2.1 Log-formulation of the Oldroyd-B system 

Let us now introduce the log-formulation proposed in [15 . We want to map solutions of the 
system psp with solutions of another system of equations where a partial differential equation 
for the logarithm of the conformation tensor is substituted to the Oldroyd-B partial differential 
equation for the conformation tensor cr. 

In order to obtain a constitutive equation in terms of '0 = In cr, following [15| , we make use of 
the following decomposition of the deformation tensor Vm G M''^'' (see Appendix [C] for a proof): 

Lemma 2.2. For any matrix Vm and any symmetric positive definite matrix cr in R'^^'', there 
exist in M.'^^'^ two antisymmetric matrices ft, N and a symmetric matrix B that commutes with 
(J , such that: 

Vu = rL + B + NcT-^. (22) 

Moreover, we have tr Vtt = ivB. 

We now proceed to the change of variable i/j — In cr. The system (fTHI) then rewrites (see [T5] 
for a proof): 

Re ( ^ +u - Vtt I = -Vj5+ (1 -e)AM+ — dive'^, 
\ot J Wi 

divi6 = 0, (23) 

^ + (it • V)-f/> ^ nt/j - + 2B + ^^{e'"^ - I)- 

It is supplied with unchanged initial and boundary conditions for u, plus the initial condition 
i/j(t = 0) = In cr(t — 0) for the log-conformation tensor. 

2.2.2 Reformulation of the free energy estimate 

A result similar to Theorem 12. II can be obtained for system (j23p . where the free energy is written 
in terms of i/) as: 

F{u, e^) = 5£ / + ^ / tr(e^ -^-I). (24) 
The following theorem then holds : 



2 ' 2Wij^ 
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e^) + (1 - e) / I Vit|2 + I tr(e^ + e"^ - 21) = 0. (25) 



V 



Theorem 2.2. Let (m,p, xp) be a smooth solution to system (j23p supplied with hom,ogeneous Dirich- 
let boundary conditions for u. The free energy satisfies: 
d 

From this estimate, we get that F{u,e'^) decreases exponentially fast in time to zero. 

Proof of Theorem (|2.2p . The proof of this theorem mimicks the proof of Theorem l2.1l We go over 
the steps of the proof, and point out the differences with the previous case. Let {u,p,ij}) be a 
smooth solution to (1^ . 

From the inner product of the momentum conservation equation in (j23|l with the velocity it, 
we obtain: 

which is equivalent to (|18p. Taking the trace of the evolution equation for the conformation tensor, 
we get: 

which is equivalent to (|2ip . Contracting the evolution equation for xjj with and using (jl3p with 
■0 = Incr, we rewrite the first term of this inner product: 



d_ 

dt 



, M-VtAj :e^^ I — +M-V)tre^ 



Recall that the decomposition (|22p of Vm allows to rewrite the second term: 

Vm : = n : + B : + (iVe^^) : e"^ = B : e'^ , (28) 

where we have used the symmetry of e'^ and the antisymmetry of ft and AT. Then, notice that, 
since and e'^ commute, we have: 

{nip - ipn) : = tr(J7i/)e^) - ti{xpne'''), 
= tr(J7i/)e^) -tr(rJV'e'^), 

= , (29) 
we finally obtain an equation equivalent to (jl9p : 

4- I tre^^2 f Vu:e^--^ f tr(e^ - /). (30) 
dtj^, U WiJt, ^ ' ^ ' 

It is noticeable that in this proof, wc made no use of the positivity of cr = e'^, in contrast to the 

proof of Theorem 12.11 

The combination dH-^WT^ (EUl+aWi^ ® gives (US): 

s [t X + 2^ X - - + " X + 2W? X + - - »■ 

(31) 

This is exactly equivalent to (fT7|) . As in the proof of Theorem 12. 11 we then obtain that F{u, e"^) 
decreases exponentially fast in time to zero. □ 
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3 Construction of numerical schemes with Scott- VogeUus 
elements for (uh^Ph) 

We would now like to build numerical integration schemes for both systems of equations psp 
and that respectively preserve the dissipation properties of Theorems 12.11 and 12.21 for discrete 
free energies similar to ([T6| and ((24|) . We first present discretizations that allow for a simple and 
complete exposition of our reasoning in order to derive discrete free energy estimates. Possible 
extensions will be discussed in the Sections 14.31 and [51 



3.1 Variational formulations of the problems 

To discretize (|15p and (j23p in space using a finite element method, we first write variational 
formulations for ([TS]) and ((25)) that are satisfied by smooth solutions of the previous systems. 
Smooth solutions {u,p, a) and {u,p, ijj) to system (fT5|) and (|23p respectively satisfy the variational 
formulations: 

= y Re + M • Vm^ • t> + (1 — e)VM : Vt> — pdivu 

+ r^cr : V'U + qdivM (32) 
Wi 

dcr \ 1 

-^+u- Vcr j : cjy - {{Vu)ct + (T{Vuf) : + — (cr - /) : 



and 



J Re (^"q^ + u ■ Vm^ ■ V + {1 — e) Vw : Vt> — p div v 



+ :^e'^ : Vv + gdivM (33) 
^ + M • VtA^ : - (OtA - ^n) : - 2B : - ^^(e-'^ - /) : 0, 

for all sufficiently regular test function {v, q, 0). 

In this variational framework, we recover the free energy estimates (fT7|) (respectively ()25p ') 
using the test functions (it,p, 2^(-f — c^^)) (respectively (^u,p, -^^^{e"^ — 1))) in ([32]) (respec- 
tively (P)) ). 

3.2 Numerical schemes with Scott- Vogelius finite elements for (uh,ph) 

Using the Galerkin discretization method, we now want to build variational numerical integration 
schemes that are based on the variational formulations (j32p and p3p using finite-dimensional 
approximations of the solution/test spaces. We will then show in the next Section 2] that solutions 
to these schemes satisfy discrete free energy estimates which are equivalent to those in Theorems 1 2. II 
and[ 
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First, the time interval [0, T) is split into Nt intervals [t", of constant size At — with 
= nAt for 71 = . . . Nt- For aU n = . . . iVy - 1, we denote by {uI,pI,(tI) (resp. {ul,pl,ij}l)), 

the value at time t„ of the discrete solutions {ufi,Ph,Ch) (resp. (ufuPhjiPh)) finite element 

spaces. 

In all the following sections, we will assume that the domain T) is polyhedral. We define a 
conformal mesh built from a tesselation of the domain V, 



made of Nk simplicial elements Kk and No nodes at the internal vertices. We denote by hx^ the 
diameter of the element Kk and assume that the mesh is uniformly regular, with maximal diameter 
h > maxi<k<NK ■ For each element Kk of the mesh 7^, we denote by tik^. the outward unitary 
normal vector to element Kk, defined on its boundary dKk- We also denote by {Ej\j = 1, . . . , Ne} 
the internal edges of the mesh 7^ when d = 2, or the faces of volume elements when d = 3 (also 
termed as "edges" for the sake of simplicity in the following). 

For the velocity-pressure field {uh,Ph), we choose the mixed finite element space (P2)'' x Vi^disc 
of Scott- Vogelius [30], where: 

• by Uh e (P2)'' we mean that Uh is a vector field with entries over D that are continuous 
polynomials of maximal degree 2, 

• and by ph £ Pi. disc we mean that ph is a scalar field with entries over 7^ that are piecewise 
continuous polynomials of maximal degree 1 (thus discontinuous over V). 

This choice is very convenient to establish the free-energy estimates at the discrete level. As 
mentioned earlier, other choices will be discussed in Section 14.31 For general meshes, this finite 
element does not satisfy the Babuska-Brezzi inf-sup condition. However, for meshes built using 
a particular process based on a first mesh of macro-elements, this mixed finite element space is 
known to satisfy the Babuska-Brezzi inf-sup condition (this is detailed in [1] for instance). The 
interest of this finite element is that the velocity field is divergence-free: 



because div Uh G Pi, disc can be used as a test function for the pressure field in the weak formulation 
of the incompressibility constraint J^divUhgh = 0. 

For the approximation of cr^ and xpi^, we use discontinuous finite elements to derive the free 
energy estimates. For simplicity, we first consider piecewise constant approximations of cr^ and xpf^ 
in Sect ions [3] and m In Section[5l we will come back to this assumption and discuss the use of higher 

d(d+l) 

order finite element spaces for CTh and i/^^. All along this work, we denote by (Th G (Pq) 2 the 
fact that the symmetric-tensor field cr^ is discretized using a '^^''^^^ -dimensional so-called stress 
field, which stands for the entries in Pq of a symmetric d x d-dimensional tensor field, thus enforcing 
the symmetry in the discretization. 




divuh{x) = 0, Va; e V, 



(34) 
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The advection terms u ■ X^cr and u ■ Vi/j will be discretized either through a characteristic 
method in the spirit of [381 or with the discontinuous Galerkin (DG) method in the spirit 

of [23] . Notice already that the characteristic method requires the velocity field to be more regular 
than the discontinuous Galerkin method in order to define the flow associated with the vector 
field Uh- 

For the discontinuous Galerkin method, we will need the following notation. Let Ej be some 
internal edge in the mesh T^. To each edge Ej, we associate a unitary orthogonal vector n = uej-, 
whose orientation will not matter in the following. Then, for a given velocity field in V that 
is well defined on the edges, for any variable in I? and any interior point x to the edge Ej, we 
respectively define the downstream and upstream values of <p by: 

(j)^{x)= lim 4>{x + 5 Uh{x)) , 

(35) 

4> {x) = lim 4>(x + 5uhix)). 

(5^0- 

We denote by |0](a;) = 4>^{x) — 4>^ [x) the jump of </> over the edge Ej and by {cj)} {x) = 
4> ix)+4' {x) mean value over the edge. Then, one can easily check the following formula for 
any function 0: 

J2 f \uH-nm^-J2[ (uh-riK,)^. (36) 



Let us now present in the next section the discrete variational formulations we will consider. 

Remark 3.1. In what follows, we do not consider the possible instabilities occurring when advec- 
tion dominates diffusion in the Navier-Stokes equation for the velocity field Uh ■ Indeed, in practice, 
one typically considers small Reynolds number flows for polymeric fluids, so that we are in a regime 
where such instabilities are not observed. 

Moreover, we also assume that < e < 1 so that there is no problem of compatibilities between 
the discretization space for the velocity and for the stress (see '61 for more details). 

3.3 Numerical schemes with cr^ piecewise constant 

Variational formulations of the discrete problem write, for all n = . . . Nt — 1, as follows: 

,n+l _n+l n+l\ 



With the characteristic method: For a given (mJJ,^;', (7^), find (it'^'+\/^'+\ cr'^'+i) e (P2)'' 
^i,disc X (Po) 2 such that, for any test function {v,q,4>) G (P2) x ^i.disc x [^o) ^ 

^ div V -\- q div ttj^^ 

ri-|-l n „ vni'-i-7i\ 



pI+^ divt; + gdiv<+i + (1 - e)Vul+^ : Vv + ^cr^' : 



At 



(37) 
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This problem is supplied with an initial condition (m^,p^,(t^) e (P2)'' x Pi, disc x (Pq) * 2^ ' . 
The function X"(t) : x E T) 1-^ X^^{t,x) E V is the "backward" flow associated with the velocity 
field mJJ and satisfies, for all x E V, 

j iX^it, x) = ul{X--(t, x)), Vt E 

\ X-{t-+\x)=x. ^^^^ 

'With. t\ie discontinuous GaZerfein method: For agiven (ujjjp^, crJJ), find (u^'^^,pJJ'^^, crJJ+^) E 
{T2f X Pi^disc X (Po)^^^ such that, for any test function q, cf)) E (P2)'' x Pi^disc x (Pq)^^^, 

O^EX^RH At 



fe=i 



' div + g div «r 1 + (1 - e) V^r^ : Vt; + ^rr^ ^ : Vt; 

c/,-((V<+iK+i+<+i(V<+if):</, (39) 



At 



+ ^ |<.n|K+il:0^ 



djd+l) 

Since cr/j G (Pq) ^ is discontinuous, we have discretized the advection term for crh with a sum 
of jumps similar to the usual upwind technique, where 0^ = (5 1^1 + {0}) (see [131 [23] ). 

Remark 3.2. In all the following, we assume that, when using the characteristic method: 

• the characteristics are exactly integrated, 

• and the integrals involving the backward flow X" are exactly computed. 

We are aware of the fact that these assumptions are strong, and that numerical instabilities may 
be induced by bad integration schemes. Hence, considering the lack for an analysis of those inte- 
gration schemes for the characteristics in the present study, our analysis of discontinuous Galerkin 
discretizations of the advection terms may seem closer to the real implementation than that of the 
discretizations using the characteristic method. 

Remark 3.3. Notice that the numerical schemes we propose are nonlinear due to the implicit 
terms corresponding to the discretization of the upper- convective derivative {Wu)cr + cr{^u)'^ 
(resp. flxp — xpfl). In practice, this nonlinear system can be solved by fixed point procedures, either 
using the values at the previous time step as an initial guess, or using a predictor obtained by 
solving another scheme where the nonlinear terms are explicited. 
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3.4 Numerical schemes with ij^i^ piecewise constant 

We now show how to discretize the variational log-formulation similarly as above. For this, we will 
need the following elementwise decomposition of the velocity gradient (see Lemma above) : 

V<+i = + Bl+^ + Nl+^e-'f'"^^' . (40) 

Moreover, for the decomposition ((40|) with u G (P2)'^, we will need the following Lemma [3.11 for 
= 1, which is proved in Appendix [Ul 

Lemma 3.1. Let ^u^^^ e i^k,disc)'^^'^ for some k € N. Then, for any symmetric positive definite 
matrix e'^h^^ £ (Pq) '' ^ \ there exist two antisymmetric matrices JIJJ^"^, TVJJ^^ e i^k,disc) ' ^ ' 
and a symmetric matrix B^'^^ G (Pfc.disc) ^~ that commutes with e^^^ , such that the matrix- 
valued function Vu^'^^ can he decomposed pointwise as: VitJ^^"^ ^h^^ + Bh^^ + N'l'^^^e^'^h'^^ . 

Variational formulations of the discrete problem write, for all n = . . . Nt — 1, as follows: 

With the characteristic method: For a given {ul,pl,'ij)l), find {ul^^ ,pl+^ ,ipl'^'^) G (P2)''x 
Pi, disc X (Po)^^^^ such that, for any test function {v,q,4>) G (P2)'' x Pi.disc x (Pq)^^^^, 

= /^Re( ^^^^p'''' +<-V<+M-i^ 



-pI+^ divt; + gdiv<+i + (1 - £)V<+i : Vv + :^e^'^^' : Vt; 



(41) 



-4(e-^^^^-/):0, 

where the initial condition (u°,p°,i/'°) G (P2)'' x Pi, disc x (Po)^^^^ is given and where X"(i) is 
again defined by 



With the discontinuous Galerkin method: For a given {u^^p^, ^^'^ (''^h^^ ' Ph^^ ' '*Ph^^) ^ 
X (Po)^^T^ such 



(Pa)'' X Pi^disc X (Po)^^^ such that, for any test function {v, q, 4>) G (P2)'' x Fi^disc x (Pq)^^^ 



= E X^R^ [ ' At + < ■ 



-pl+Uivv + qdivu]'+^ + {l-e)Vul+^ ■.Vv + ^^e'''"^*' : 

At j'-"^' ("r vr^ - ^r^^r^) : - ■. (42) 



Ne „ 

/ |<.n|I^ri]:c/>^ 

..-1 Je. 



= 1 -"^3 
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3.5 Local existence and uniqueness of the discrete solutions 

Before we show how to recover free energy estiraates at the discrete level, let us now deal with the 
local-in-time existence and uniqueness of solutions to the discrete problems presented above. 

First, since the mixed finite element space of Scott- Vogelius chosen in the systems above for the 
velocity-pressure field satisfies the Babuska-Brezzi inf-sup condition, notice that the system ([57]) 
is equivalent to the following for all n = . . . Nt — 1: For a given (mJJ;, cr^), find {u^~^^, o"^^^) G 
(IP2)lv=o X in)^ such that, for any test function {v,cj)) G (P2)L=o x i^o)^ , 

"-^ : - ((V«r^)<+i + al^Wut^f) : c/> (43) 

1 



where the flow is defined by ((55)1 and where we have used the following notation: 

(P2)div=o = e J^qdWv^O, ^qe Pi, 



,disc 



(44) 



Notice that it is also straightforward to rewrite the systems ([5^ . ([1T|) and ([1^ using G 
(P2)div=o instead of [uh^Ph) & (^2)'^ x Pi, disc- For instance, the system (|1T|) is equivalent to: For 
a given K,^^), find G (P2)l,=o x i^)^ such that, for aU (t;,0) G (P2)^i,=o x 



, d(d+l) 

Re ( ^ + < • Vul+^ j-v+il- e)Vul+^ : Vv + ^e^^ : Vt; 

-^(e-<"-/):0. 
Then, we have the: 

Proposition 3.1. For any couple (mJJjCtJJ) with af^ symmetric positive definite, there exists cq = 
Co (M5J,crJJ) > such that, for all < At < cq, i/iere exists a unique solution (mJJ^^, crj^^^) to t/ie 
system ([37]) (resp. p9p ) to^/i crJJ^^ symmetric positive definite. 

Proof of Provosition \3.1[ The proofs for systems ([57]) and ([55)1 are fully similar, so we will proceed 
with the proof for system ([57]) only, using its rewriting as system (|15|) . 

For a given mesh 7/j, let us denote by g M.'^^o+3Nk ^^j-^g vector whose entries are respec- 

tively the nodal and elementwise values of (mJJ^^, crj^^^), solution to the system The system 

of equations ((431) rewrites in terms of the vector y"+^ G M.'^^o+sNk g^. Jq^. given y" and At, 
find a zero 1""+^ of the application Q defined by 

Q{At, y"+i) = AM(y"+i)y"+i + AiB(y")y"+i y"+i - c(y", At) , (46) 
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where A and B are linear continuous matrix-valued lunctions in ]]j(2JV£j+3A'ir)x(2Af£,+3WA')^ a^fj 
where C is a vector- valued lunction in ^^Nd+^Nk (notice that the dependence of the function C 
on At is only related to the computation of the backward flow during a time step At, so that 
C(r",0) = F", and with the DG method it simplifies as C(y", At) = y"). The functions A, B 
and C also implicitly depend on 7^, as well as on the parameters Re, Wi, e. 

Now, Q(At, Y) is continuously differentiable with respect to (At, Y) and we have, with / the 
identity matrix in m(2Wi.+3Wk)x(2W^+3Wk). 

Vi'Q(At, Y) = I + At{B{Y'') + A{Y) + (Vy A)y) . (47) 

Then, for given vectors F" and Y, the matrix Vy(5(At, Y) is invertible for all < At < 
\\B{Y'') + A{Y) + (Vyyl)y|P^ (with convention \\B{Y'') + A{Y) + {VyA)Y\\^^ = oo if B(y") + 
yl(y) + {WyA)Y = 0), and then defines an isomorphism in ]r2Wd+3JVk^ 

Let us denote by S'^ the subset of r'^^d+sNk f^i^.^^^ Q^ly contains vectors corresponding to 
elementwise values of positive definite matrix-valued functions (Th in D. Since 5^(R''^'^) is an 
open (convex) domain of M''^'^, S^^ is clearly an open (convex) domain of . 

Since Q{0,Y") = and Vy(5(0,y") is invertible, in virtue of the implicit function theorem, 
there exist a neighborhood [0, cq) x V{Y"-) of (0, y") in M+ n S'^ and a continuously differentiable 
function R : [0, cq) -> V(Y"), such that, for all < At < cq: 

Y = R{At) <=> Q(At, y) = 0. 

For a given time step At e [0, cq) and a given symmetric positive definite tensor field crJJ, i?(At) e 
V{Y^) is the vector of values for a symmetric positive definite solution o"J^^^ to the system pS)) . 
Notice that, up to this point, cq — co(y") is function of y", as well as Re, Wi, e and 7/i. □ 

For solutions {u}l,(7fj to the systems (j4T|) and (|42l) . we similarly have: 

Proposition 3.2. For any couple {u"^,xp]t^), there exists a constant cq = CQ^u^jip^) > such that, 
for all < At < cq, there exists a unique solution {u^'^^ to the system (|4ip (resp. (|42p j. 

The proof of Proposition l3 . 2l is fully similar to that of the Proposition l3.11 but for the expressions 
of Q(At,y) with respect to Y. An additional term AtZ?(y) now appears in Q due to e^^* . 
This term is continuously differentiable with respect to y, and the derivative Vy'(5(0,y") is still 
invertible. Thus, the proof can be completed using similar arguments. 

Anticipating the results of Section [HI we would like to mention that the above results will be 
extended in two directions, using the discrete free energy estimates which will be proved in the 
following. 

• We will show that the constant cq in Proposition 13. II (resp. Proposition 13. 2p can be chosen 
independently of (itJJ, crJJ) (resp. (mJJ, iph))^ which yields a longtime existence and uniqueness 
result for the solutions to the discrete problems (see Propositions 16.11 and 16.21 below) . 

• We will also show, but for the log-formulation only, that it is possible to prove a long-time 
existence result without any restriction on the time step At (see Proposition 16.31 below) . 
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4 Discrete free energy estimates with piecewise constant 
discretization of the stress fields CTh and 

In this section, we prove that various numerical schemes with piecewise constant (Th or xp/^ satisfy 
a discrete free energy estimate. We first concentrate on Scott- VogeHus finite element spaces for 
{uh,Ph) and then address the case of other mixed finite element spaces in Section [ 



4.1 Free energy estimates with piecewise constant discretization of cTh 
4.1.1 The characteristic method 

Proposition 4.1. Let (itJJjpJJ, cr)J)o<„<7VT be solution to p7p . such that crJJ is positive definite. 
Then, the free energy of the solution (u^,p^,(T^): 



(48) 



satisfies: 



F. 



n+1 



Re. 



ul+^ - u 



n\2 



V 



+ At[ (l-e)|VM 
Jv 



n+l|2 
h I ■ 



tr«+i + (^ri)-i-2/)<0. 



(49) 



In particular, the sequence {F^)Q<n<NT ^•s non-increasing. 

Proof of Proposition]^ Let (u^+\p^+\ cr'^'+^) be a solution to system Notice that (cr'^'+^)-i 
is still in (Po)^^^^^. We can thus use (itJJ+^pJJ^^ (/ - {(t]['^^)~^)) as a test function in the 
system (|57| . which yields: 



Re 



I n+l|2 



Pl+' div<+i div<+i + (1 - e)|V<+Y 



^n+l 



Wi 



n+1 



2Wi 



At 



(50) 



{I -{err')'') 
1 



2(V<+i)<+i : (J - «'+i)-i) + — (<t;:+i - /) : [I - {^1+^)-^] 



We first examine the terms associated with momentum conservation and incompressibility. We 
recall that wJJ^^ satisfies ([M)) since we use Scott- Vogelius finite elements. By the Stokes theorem 
(using the no-slip boundary condition), we immediately obtain: 



/ <.vK+ip = - / (div<)K+ip = o. 

Jv Jv 
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The terms involving pJJ^^ also directly cancel. We now consider the terms involving The 
upper-convective term in the tensor derivative rewrites: 

(V<+i)<Tri : (7 - = : V<+i - div<+\ 

which vanishes after combination with the extra-stress term cr^^^ : Vm^^'^^ in the momentum 
conservation equation, and using the incompressibility property. The last term rewrites: 

K+i -!):(/- K+^)-^) - tr «+i + - 21) . 

The remaining term writes: 

/ K+i-<T5:oX"(r)):(J-K+i)-i)= / trK+i)-trKoX"(r)) 

+ tr(KoX"(t«)][a;j+i]-i-/). 
We first make use of ^ with (t = (tI o and r = (tI+^: 

tr (K ° X"(r)]K+i]-i - J) > trlnK ° ^"(^")) - trln«+i). 

Then, we have: 

/ - tr [al o + Hal o = / - tr (<t^ + Inrr'^^) , 

Jv Jv 

since the strong incompressibility property (divitJJ — 0) implies that the flow defines a 

mapping with constant Jacobian equal to 1 for all t £ [i", Finally, we get the following lower 

bound: 

/ K+i-^5JoX"(r)):(J-K+i)-i)> / trK+i-ln<+i)-trK-lnO, 
Jv Jv 

hence the result (P^ . 

Notice that tr (cr^+^ + (cr,"+^)"i - 21) > in virtue of the equation which shows that the 
sequence (-F)")o<n<A'T is non-increasing. □ 

4.1.2 The discontinuous Galerkin method 

Proposition 4.2. Lei (mJJ,^^^, (TjJ)o<„<Ary &e solution to p9p . smc/i t/iat crJJ zs positive definite. 
Then, the free energy defined by (|48p satisfies the free energy estimate ()49p . /n particular, the 
sequence {FJ^)o<n<NT non-increasing. 

Proof of Proposition \4.^ We only point out the differences with the proof of Proposition l4.1l They 
consist in the treatment of the discretization of the advection terms for cr^. We recall that the test 

Ne 



function in stress is (p — 2^iil ^ i^h^^) so that we have 



Y^f \ul-n\lal+^:{l- 

Ne „ Ne ^ 

|<.n|ItrK+i)l+^ / |<.n|tr(^ri-K+i-+)-i-j). 



j=l " j=i ■ 



INRIA 



Dissipative Oldroyd-B schemes 



21 



Again, we make use of PU)) . with cr = cr^^^' and r = crj^^^'^: 

We get, by Formula ([55]) . the fact that e (Pq) ' ' , the Stokes theorem and the incompress- 

ibility property ((M)) : 



^ / K • n| K+i] : (7- K+i)-i)+ > ^ / \ul ■ n\ [trK+i - Ina^ 

Nk „ 

= -5] / K-n^JtrK+i-ln^+i), 

Nk „ 

= -5^(trK+i-h,<+i)) / 

k=l -^9^" 
Nk „ 

= -5](trK+i-ln<+i)) / divK), 

= 0. (51) 

Moreover, it is easy to prove the following, using the same technique as in the proof of Proposi- 
tion O 

/ K+i-^','):(/-K+i)-i)> / trK+i-ln^ri)-trK -lnO. 
jt> Jt> 

This concludes the proof. □ 

4.2 Free energy estimates with piecewise constant discretization of -0^ 

This section is the equivalent of the previous section for the log-formulation. 

4.2.1 The characteristic method 

Proposition 4.3. Lei (ttJJ,^^^, ■05J)o<n<Arr &e solution to (PTj ). Then, the free energy of the solution 
= FK,e«) ^ Kf + ^J^tr(e« - - I) , (52) 



satisfies: 



At ^^(l-e)|V«r^|^-|-^tr(e^;:^^-fe-^;.'^^-2/) <0. 



(53) 



/n particular, the sequence {F^)o<n<NT ^•s non-increasing. 
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Proof of Proposition\4^ We shall use as test functions (ul'^^,pl+^, IWi (^'^''^' ~ (SB)- 
emphasize that, as long as the solution (itj^^^jp^'^^, -05^^^) exists (see Proposition 13. 2p . e^^^ is 
well-defined, symmetric positive definite and piecewise constant. 

The terms are treated similarly as in the proof of Proposition 14.11 For the material derivative 
of xp,^, we have: 



+1 



•D 



Using the equation (fTT|) with a = e'^'>^ and t — e'^'^°^ \ we obtain: 

-iplo : e'^^^' > tr(e^^^' - e^^^"(*")), 

and thus: 

(tA^i - xPl o : (e^^' - J) > / trle'^^^' - V^') - / tr(e^^' - xj^l) o 



where the fact that the Jacobian of the flow X" is constant equal to one (because is divergence 
free) has been used in the change of variable in the last equality. 
Besides, using the equation (P^ . we have: 

Jv Jv 

= 0. 

Last, using f^ : 

Jv Jv Jv 

= f Vul+' : e^^^' - / div«+i), 
Jv Jv 

= I :e^^^\ 
Jv 

which cancels out with the same term J.^ e^^'^ : VitJ*^^ in the momentum equation. □ 
4.2.2 The discontinuous Galerkin method 

Proposition 4.4. Let (mJJ,pJJ, ■0JJ)o<n<JVr solution to (H^ . Then, the free energy FJ^ defined 
by ([5^ satisfies the free energy estimate ([55]) . In particular, the sequence (-F)")o<n<JVT non- 
increasing. 

The proof is straightforward using elements of the proofs of Proposition [¥31 and Proposition^ 
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4.3 Other finite elements for {uh,Ph) 

In this section, we review some finite element spaces for {uh,Ph) other than Scott- Vogelius for 
which the results of the last two sections still hold. 

First, let us stress the key arguments we used in the proofs above. If the advection terms 
u ■ Vcr and u ■ V'0 are discrctizcd by the characteristic method, we need the velocity field itJJ to 
be divergence free: 

div< = 0, (54) 

in order for the flow X" satisfying ([38|l to be with Jacobian one. When is only piecewise 
smooth (consider below the case of Pi, disc velocity fields), the divergence in the left-hand side 
of ([5^ should be understood in the distribution sense. By the way, the equation univoquely 
defines the trace of the normal component itJJ • n on the edges of the mesh, which is a sufficient 
condition to define the flow associated with the vector field mJJ through ([55]) . and which is necessary 
to treat the advection term in the Navier-Stokes equation. 

If the advection terms are discretized by the discontinuous Galerkin method, it is necessary 
that the trace of the normal component of Uh be univoquely defined on the edges of the mesh since 

it appears in the jump terms Ef=i Ie, \K ' "I K^'l ■ Ef=i Ie, \K ' "I 1'/''^'! ■ 

the variational formulations. But to obtain (|5ip . and contrary to the characteristic method, only 

the following weak incompressibility property is needed: 

\fk=l...NK, I div< = 0, 

which is equivalent to write 

yq e Po, / div«)g = 0. (55) 
Jv 

These properties needed to obtain the discrete free energy estimates are summarized in Table [T] 



Advection discretized by: 


Characteristics 


DG 


Requirements for Uh'- 


div Uh — 

( ^ det(V,X") = 1 ) 

( {uh ■ n)\Ej well defined ) 


f^qdivuh = 0, Vg e Po 
and 

{uh ■ n) \Ej well defined 



Table 1: Summary of the arguments with {uh^Ph, fh) or {ufi,Ph, '4'h) (P2)^' x Pi. disc x (Po) ^ 



Below, we consider the following alternative choices of the finite elements space for {uh,Ph)- 

• the Taylor-Hood finite element space: {uh,Ph) G (P2)'' x Pi, which satisfies the Babuska- 
Brezzi inf-sup condition, whatever the mesh ; 

• the mixed Crouzeix-Raviart finite element space (see |12,): {uh,Ph) G {Pi^)'^ x Pq, where 
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which also satisfies the Babuska-Brezzi inf-sup condition, whatever the mesh ; 
• stabihzed formulations for {uh,Ph) G (Pi)'' x Pi or {uh,Ph) G (Pi)'' x Pq. 

This is not exhaustive, but it is sufficient to highlight which modifications are needed in the 
variational formulations, compared to the Scott- Vogelius mixed finite element, for the discrete free 
energy estimates to hold. In particular, some projection of the velocity field is needed in the 
discretization of the advection terms u ■ ^cr and u ■ Vt/^ in order to satisfy the requirements of 
Table [H These projection operators are introduced in the next Section 14.3.11 
The results of Section 14.31 are summarized in Table [2l 

4.3.1 Some useful projection operators for the velocity field 

Let us introduce three projection operators for the velocity field. 

Following [3H], we first define the orthogonal projection P™* of (Pi. disc)'' onto the piecewise 
constant solenoidal vector fields built from affine continuous scalar fields: 



We suppose here that c? = 3, but the extension to the case d = 2 is straightforward. We set 
P™* (m/i) = V X where V'/i € (Pi)'', such that -0/1 x n\Qx> = 0, satisfies : 



Because P™* (m/i) is solenoidal, we always have the strong incompressibility property ([54|) : 



for any velocity field u^. Of course, this operator is consistent only for divergence free velocity fields 
U}i (or velocity field tt;,, with vanishing divergence when h goes to zero). See for consistency 
results. 

Second, following [13], we define the Raviart-Thomas interpolator P^^° from (Pi, disc)'' onto 
the vector subspace of (Pi. disc)'' made of the vector fields in (Po)'' + JcPq with continuous normal 
component across the edges Ej (whose trace on Ej is then univoquely defined). The projection 
P^'^° (mJJ) clearly satisfies, for any element Kk'- 



{V X 0.|a e (Pi)'', 0> X n = on dV}. 




divP™* [uh)=Q 




k 1 



(57) 



Thus, it satisfies the weak incompressibility property fSS]) : 
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if, and only if, the velocity field uf^ also satisfies it. 

Third, we define p^^'^^ as the Brczzi-Douglas-Marini projection operator [71[H|- It is with value 
in (Pi)''. This projection operator satisfies the same divergence preservation property ((57)) than 
P^'^° , but is of better accuracy. 

Note that p^^^^ like P^'^° are local interpolating operators in the sense that all the computa- 
tions can be made elementwise. This is not the case for P™*. 

In addition, we will need the following lemma: 

Lemma 4.1. For any velocity field such that the previously defined interpolating operators are 
well defined, the normal components of the interpolated vector field, P[°* (u^) ' n, P^"^" (u^) ■ n 
and pBDM (njj) . fi (jj-g algo well defined on any internal edges Ej. 

Moreoever, if u^,^ G (Pi, disc)'' is a velocity field such that, for all k ^ 1 . . . Nk'- 



then div{P^'^° (uD) = div(P,f (m^)) = (in the sense of distribution). 
Proof By construction, P^'^" and pBDM 

are with value in velocity fields such that their normal 
component is continuous across the edges. This is also the case for P™* since P™* is with value 
in divergence free velocity fields. Then, from the equation (fST)) . we have div(P^"^° (ufj) = 0. 
Since div(P^'^° (u^J)) is in {Fof, this shows that div(P^'^« (mJJ)) is zero in any element Kk ■ Finally, 
Pj^^° {ufj has continuous normal components across the edges of the mesh. This shows that 
div(P^"^'' (u^j^)) = in the sense of distribution. The same proof holds for the projection operator 



4.3.2 Alternative mixed finite element space for {uh,Ph) with inf-sup condition 

In this section, we show how to derive discrete free energy estimates with mixed finite element 
spaces for the velocity and pressure fields which satisfy the inf-sup condition, but which are not 
the Scott- Vogelius finite elements. 

Let us first consider the Taylor-Hood element for {uh,Ph), that is (P2)'' x Pi. In this case, since 
the velocity field Uh is not divergence free either in the weak form (|55p , or in the strong form (|54p , 
a projection of the velocity field is required in the discretization of the advection terms u ■ 
and u ■ Vi/j. More precisely, we need to use the projection velocity PfJ'^u^ (and, among the three 
projection operators we introduced above, this is the only one which is such that the strong or weak 
incompressibility is satisfied). For the characteristic method, one uses the flow X"'{t) satisfying: 




>BDM 



□ 



ix^{t,x) = P-°'ui{x^{t,x)), e 

X"(t"+i,a;) = X. 



(58) 
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For the discontinuous Galerkin method, the advection term in the conformation-tensor formulations 
writes (see the last line in (15^ ): 

Ne „ 

/ iprK)-r^iK+^i:0+. 

Notice that in the terms |o"^"''^] : 0^, the projected velocity P^^u)^ is used to define the upstream 
and downstream values following ([55)1 . Another modification, which is specific to the Navier-Stokes 
equation, is needed to treat the advection term on the velocity. Namely, one needs to add to the 
weak formulation the so-called Temam correction term (see j41)): 

+ — JjW{K){vul+') (59) 
in such a way that, when u'^^^ is used as a test function: 



Re / {K ■ V<+i) • <+i + ^ / div (K) \u 



With these modifications (projection of the velocity field in the advection terms, and Temam cor- 
rection terms), the free energy estimate (|^^ is satisfied by the scheme. Similar results (discrete free 
energy estimates for {ufi,ph,ipii) in (P2)'^ x Pi x (Pg) ' ' ) can be proved on the log-formulation. 

Let us now discuss the use of Crouzeix-Raviart finite elements for velocity: (uh,PhT(^h) in 
gpCR^jd X Pg X (Po)^-^ (see In this case, the Navier-Stokes equations can be discretized 

using a characteristic method: 



, 1 J Is 



Re 



- divv + q div + (1 - £) Vm"+^ : V^; + -^f 



where X" is obtained from the projected velocity field P/imJJ as: 

f ^X"(t) - P^<(X"(i)), yt e [t",t"+i] 

I X"(t"+i) = a;. 



(61) 



The projected velocity Phu'h defined using any of the three projectors presented above, that 
is Pj^^^u'f^ or P^^^^u^. The Navier-Stokes equations can also be discretized using a 

discontinuous Galerkin formulation: 

0=Ey^ Re ( " ' + «) ■ ^ J ■ + Re ^ |P,«).n|K+i].M 



fc=i 
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Here again, PhU^ is any of the three projectors presented above. We would hke to mention that 
we are not aware that the projector P™* has ever been used with discontinuous Galerkin methods, 
so that the consistency of the discontinuous Galerkin approach combined with this projector still 
needs to be investigated. 

Likewise, the advection term u ■ in the equation on the stress can be treated by the 
characteristic method or the discontinuous Galerkin method, as above for the advection term in 
the Navier-Stokes equations. 

Notice that whatever the projecting operator used, drv^PhU^) = holds (see Lemma ini above). 
With this property, it is easy to check that Propositions l4. 1 1 and 14. 21 still hold for this finite element. 
For example, the advection term in the Navier-Stokes equations is treated as follows (using the 
fact that dW{Ph (u^^)) = and ([551): 




= 0. 



Discrete free energy estimates for {ufi,Ph, tph) (Pf "'^■)'' x Pq x (Pg) ' 2 ' can be similarly proven 
on the log-formulation. 

4.3.3 Alternative mixed finite element space for {uh,ph) without inf-sup 

It is also possible to choose a mixed finite elements space for {uh,Ph) that does not satisfy the 
Babuska-Brezzi inf-sup condition, like (Pi)'^ x Pq or (Pi)'' x Pi. The loss of stability due to the 
spaces' incompatibility can then be alleviated through a stabilization procedure, like Streamline 
Upwind Petrov Galerkin, Galerkin Least Square or Subgrid Scale Method (see [3TJ [TTl (TH] ) . In the 
following, we consider very simple stabilization procedures, for which only one simple quadratic 
term is added to the variational finite element formulation in order to restore stability of the 
discrete numerical scheme. 
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Let us first consider the mixed finite element space (Pi)'' x Pq for {uh,Ph)- If the term u ■ 'Vcr 
is discretized with the characteristic method, the system then writes: 



Nk 



At ■ n ' ,1 J '2 

-pI+^ dWv + qdivul+^ + (1 - e)V<+i : Vt; + ^^cr^^ : Vt; 



"i;"""^'"^ ) : </> - ((V«r^) + (V<+^)") : (63) 



1 

^ Wi 

We 



with a flow X" computed with the projected field P™* (itJJ) through ([58|) . The projection operator 
is the only one we can use among the three projectors we introduced in Section [4.3.11 because 
the weak incompressibility property (|55p is not satisfied by mJJ. 

The stabilization procedure used in has been studied in [3S]. Proposition 14.11 holds for sys- 
tem ()63p . its proof being fully similar to the case of Taylor-Hood finite element (see Section l4.3.2p . 
since the additional term Xlj!fi \ /e IP'jIM is non negative with the test function used in the 
proof. All this also holds mutatis mutandis for discretization of the advection terms by a discon- 
tinuous Galerkin method, and for the log-formulation. 

Let us finally consider the mixed finite elements space (Pi)"* x Pi for {uh,Ph)- If the term u 'Vcr 
is discretized with the characteristic method, the system then writes: 

= / Re f < + ul ■ VulA ■v + ^ div<(^ . <+i) 



V 



At " " J 2 

-p'^'+i div-u + qdivu^'+i + (1 - £)V<+i : Vv + 4f<^^ ■ '^^ 



(TloX"{t") 



At 

Wi^ ' ^ 



Wi 

n+l\ 



k=l 



Nk 

with a flow X" again computed with the projected field P^°* (itJJ) through ([58|) . Again, we are 
led to choose the projection operator P™* because the weak incompressibility property ([55]) is not 
satisfied by m)^. The stabilization procedure used in has been studied in [S]. Proposition 14. II 
holds for system , its proof being fully similar to the case of Taylor- Hood finite element (see 
Section [4.3.2^ . since the additional term J2k=i ^k^ ^Ph ' ^1 ^^'^^ negative with the test 
function used in the proof. All this also holds mutadis mutandis for discretization of the advection 
terms by a discontinuous Galerkin method, and for the log-formulation. 
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• Advection discretized by: 


Characteristics or DG. 


• (uh , Ph) in . . . , 


=^ equations modified: 


Srntt-Vop'pliim 


( Tint n 1 n c 1 






T'a vl or-TToorl 


1 prot 


(P2)'' X Pi 


+ Temam term 


Crouzeix-Raviart 


_^ pBDM pRTo prot 




+ Ph (mJJ) for Navier term 


Stabilized (Pi)'' x Pi 


1 TDrot 

+ Temam term 


stabilized (Pi)'' x Po 


1 proi 

+ Temam term 



Table 2: Summary of some possible finite elements for {uh,Ph,o'h/tPh) when a-h/il^h ^ (^^o) ^ 
with some possible projections for the velocity field (see Section . 



5 Higher order discretization of the stress fields cr/j and 

We now show how to build numerical schemes with higher order discretization spaces for the stress 
that still satisfy a discrete free energy estimate. We typically have in mind piecewise linear spaces 
for cT/j and "0,^. 

From the previous proofs establishing discrete free energy estimates, it is clear that we need 
to use nonlinear functional of cTh and ■0;^ as test functions, namely cr^^ and e^f^. Finite element 
spaces other than Pq are typically not invariant under such nonlinear functionals, and this brings 
us to introduce projections of these nonlinear terms on Pq. A Po-Lagrange interpolation operator 
T:h is convenient, because it commutes with nonlinear functionals (see Lemma 15.21 below). 

Choosing a Po-Lagrange interpolation operator for test functions will also have another impor- 
tant consequence, because it will require that Pq be a subspace of the finite element space used to 
discretize the stress. And our proofs establishing discrete free energy estimates will use in a crucial 
way the fact that the interpolation operator coincides with a orthogonal projection onto Pq (see 
Lemma TS . 1 1 below) . The need for iTh to coincide with a orthogonal projection onto Pq also a 
priori limits the maximum regularity of the discretization of the stress, essentially to piecewise Pi 
finite elements. Therefore, we consider CTh and ipy^ in either of the following finite element spaceil|: 

^ <i(ti+l) , ^ d(d+l) 

(Pl+Po)^ or {¥l.d^sc)^ ■ 

In Section IS.l) we introduce the interpolation operator tt/j. Then we prove that, for a Scott 
Vogelius discretization of the velocity-pressure field, a free energy estimate can be obtained for 

d(d+l) d(d+l) 

■'Note that, clearly, (Pi -|- Pq) 2 is only a subspace of (f'x^diac) ^ 
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discretization schemes close to those considered in Section 21 when ah (respectively ■0^) is in 
(Pq) * 2* ' . This is the purpose of the Section [521 (respectively Section 15. 3p . Finally, we show 
in Section 15.41 how these results can be extended to other finite element discretizations of the 
velocity-pressure field. 

5.1 The interpolation operator vr^ 

Let us introduce the projection operator tt/,, as the Pq Lagrange interpolation at barycenter Ok^ 
for each Kk 

Definition 5.1. For k — 1. . .Nk, we denote by Ok^ ^he barycenter of the triangle K^. For any 
(p such that Vfc — 1...Nk, <P{(^Kk) well-defined (for example (j) is a tensor-valued function, 
continuous at points Ox^j^ we define its piecewise constant interpolation by : 

yk = \...NK ,Tlh{(f>)\K,=(t>{OK,)- 

Notice that this definition also makes sense for the case in which is matrix- valued. And 
this interpolation operator Hh coincides with the L orthogonal projection from (Pi, disc) ^ onto 

d(d+l) , , ... I I 

(Pq) 2 , as shown in the following Lemma proved in Appendix [Dl 

Lemma 5.1. Let tt^ be the interpolation operator introduced in Definition \ 5.1[ Then, for any 

d(d+l) 

4>h e (Pi, disc) 2 , we have : 

f ct^h ■■ = I {(f>h) ■■ 4, V4 e (Po)"^. 
Jv Jv 

In addition, the following property holds, which is important in the choice of this particular 
interpolation: 

Lemma 5.2. Let tt^ be the interpolation operator introduced in Definition \5.1\ The interpolation 
operator vr/i commutes with any function f: for any functions f and cf)^ such that (p^ o-nd f{4>h) 
are well-defined at the bary centers 6k, 

(/(0J) = /(tt,, (0J). 

The proof of Lemma 15.21 is straightforward since, by Definition 15.1) the interpolation iTh only 
uses specific values at fixed points in the spatial domain T>. 

5.2 Free energy estimates with discontinuous piecewise hnear cTh 

In this section, we consider the following finite element discretization: Scott- Vogelius (P2)'' x Pi, disc 
for (uh,Ph) and (Pi,disc)^^^^ or (Pi -FPq)^^^^ for cru- 
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5.2.1 The characteristic method 

If the advection term u ■ ^cr is discretized by the characteristic method, the system writes: 
- pI+' divv + gdiv<+i + (1 - £)V<+1 : Vv + ^^^^ {ctI+^) : 



At 



Wi^ /I > ^' 



(65) 

where X" is defined as in Notice that we have used the projection operator tt^ in four terms. 
It will become clearer from the proof of the free energy estimate below why those projections are 
needed. 

Proposition 5.1. Let crJJ)o<„<jVr he solution to (|65p . such thatiTh (""/J) is positive definite. 

Then, the free energy of the solution (u^,p^,<T^): 

FJ: = FK, TT, K)) ^ ^ + _1_ ^ tr K [al) - Iutt, «) - I) , (66) 



V 



At |^(l_e)|V<+i|2 + ^tr(^,«+i)+7r, «+i)-'-2/) < 



(67) 



In particular, the sequence (i^^')o<n<Afr '•s non-increasing. 

Remark 5.1. The ensemble of symmetric positive definite matrices is convex. This implies that 
a piecewise linear tensor field is symmetric positive definite as soon as it is symmetric positive 
definite at the nodes of the mesh. Moreover, this also implies that iTh (cTh) is symmetric positive 
definite as soon as cr/j is a piecewise linear (possibly discontinuous) symmetric positive definite 
tensor field. 

Remark 5.2. It is easy to extend the result of Proposition CO] to show that, if At > is small 
enough and crJJ is positive definite, then there exists a unique solution to with cr^^^ positive 
definite. 



Proof of Froposition \5.1\ The test functions we choose are ^ujj^^jpj^^^, 2^(-f — tt/i (crJJ^^) 

Recall that by Lemma f5.2| (tt^ (ctJ^^^)) ^ = tt/i ((cJ^^^)^^) • The proof is similar to the one of 
Proposition 14. II except in the treatment of the constitutive equation. 



RR n° 6413 



32 



Boyaval et al. 



The upper-convective term in the tensor derivative writes (using Lemma 15.11 and the incom- 
pressibihty property ([M)) ): 




which vanishes after combination with the extra-stress term in the momentum equation. 
The last term rewrites (using again Lemma 15.1^ : 




The remaining term writes (using Lemma fS.ll the equation (fTU|) with a = tt/, (<t}J) o X"(i") 
and T = TTh (""h"''^); and the fact that the Jacobian of X" remains equal to one due to the 
incompressibility property (134^ '): 

/ {al+'-7:n{cTl)oX-{e)) : {I - i^n {<tI+')~') 
Jv 

= - iTTTu {ctD o X"(r) + tr (vr,, (al) o X^{e)^h - /) , 

> / tr<+i - tr^^ {ctD o X"(r) + trhiTT,, (^,'J) o X"(r) - trhw,, , 

JV 

= [ ir^n{cTl+^)~tr^h{cTl)+ir\niih{rTl)-iT\n^h{cTl+^). 
Jv 

This completes the proof. □ 
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5.2.2 The discontinuous Galerkin method 

If the advection term u ■ ^cr is discretized by the discontinuous Galerkin method, the system 
writes: 



k=l 



At 

■ p]l+^ divt; + gdiv w^'+i + (1 - e)V<+i : Vv + ^^tth K+^j : Vv 



At 

Ne „ 

/ K.n|[^,K+i)]:0^ 



As for the characteristic method, the projection operator tt/i is used in four terms. Besides, hke in 

d(d+l) 

the case where crh g (Pq) 2 ^ the advection term u ■ Vcr is discretized using a jump term only. 
Indeed, in order to derive discrete free energy estimates, we treat the discrete advection term using 
the projection Wh [crh) G (Pq) ^ of the stress field cr^, the derivative of which is zero. 

Proposition 15.11 still holds for the system ([55]) . The proof is straightforward using all the 
arguments of the previous sections, except for the treatment of the discrete advection term for 
u ■ Vcr. Using the equations (fT(I|) . ((5^ . the fact that tt/j (cJJ^^) g (Pq) * ' and the weak 
incompressibility property ([55|) . we have: 

Ne 



5^ / K . n| K+i)] (j - TT, 

j=l 

Ne j. _i 

= / K.n|[tr7r4<+i)] + K.n|tr(7r,(^;:+i--)7r,(^ri-+) - /) 
>^ / K.n|[tr7r4<+i)] + K.n|tr(ln^,(,T;:+i'-)-ln7r,, (^;:+i 

= J2 / K ■ "I Mtt, - luTT, K+i))i, 
fc=l -^^^^ ^ ^ 

= 5] - tr (tt, (^^1) - luTT, {crl+') ) / divK) , 



fc=l 
=0. 
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5.3 Free energy estimates with discontinuous piecewise linear ij^i^ 

In the following section, we write free-energy-dissipative schemes using the log-formulation with i}}f^ 
piecewise linear. For this, we again need the projection operator Tih introduced in Definition 15.11 
We consider the Scott- Vogelius finite element space for {uh,Ph) and the following decomposition 
of the velocity gradient Vit/j e (Pi, disc) ' ' : 

Vuh = n,, + Bh + NhTTh {e^-y' . (69) 

Notice that since (e'^'^)"^ g-'^'-f'^'^) is in (Pq)^^^, we have Vlh,Bh,Nh e (Pi.d^^c)^^^ in 
virtue of Lemma 13.11 with fc = 1 . 

5.3.1 The characteristic method 

If the advection term u ■ ^cr is discretized by the characteristic method, the system writes: 
= / Re ( <^\-< + ul ■ Vm"+i 



pI+^ div-u + gdiv<+i + (1 - e)V<+i : Vv + — tt,, [e^"^ j : Vv 



At 



(70) 



In the system above, we have used the projection operator tt^ to treat the same terms as in the 
system (p5)) . But in addition to these, we have also used the projection operator for the exponential 
term e~'^'>^ in the Oldroyd-B equation. 

Proposition 5.2. Let {u^,p^,xp]^)Q<cn<NT &e solution to (|70p . Then, the free energy of the solution 



Fl}^F(u^t.e--»'^^]^^ I + / tr 



2 ./p' ' 2Mji, 



, (71) 



satisfies: 



(72) 



+ Mj (l-e)|V<+Y-f-^^tr(e'^''(^J^)+e-'^'-(^")-2l) < 0. 
In particular, the sequence {FJ^)(Xn<NT '■^ non-increasing. 

Proof of Provosition 1 5. Si The proof is similar to that of Proposition l4. 31 except for the terms using 
the interpolation operator tt^. We shall use as test functions 2Wi(^'» ^e''"'"^ 

in (|^T|) . And we will make use of the following property all along the proof (see Lemma E 

nn{e^-^')=e^-i^-^'). 



INRIA 



Dissipative Oldroyd-B schemes 



35 



For the material derivative of tpf^, using Lemma 15.11 the equation (llip with cr — e'^^ and 
T — e^ft°-'^"(*"), and the fact that the Jacobian of the flow X" is one for divergence free velocity 
field mJJ, we have: 

^ - TT, {^D o X"(t")) : (vr, (fl'') - /) 

= / (vr. (^r^) - (^;:) o : e--^^l^') tr (tt, (^^i) - (t^^J) o , 

J-D 

> ^ tr (e'^^^r - TT,, - ^ tr (e^^^^^) - (t/.',')) o 

= ^tr _ (t^r^)) - ^tr [e^-^^-'> - (.^D) ■ 

Besides, using the equation ([^ . we have: 

/ (jirv, (v^r^) - (t^r^) f^r^) : ^ej^-^^i^^) - 1) 

J-D 

= 0, 

and using the equations and (|34p : 

: (tt, (e^^'"^) - /) = / : e^-i^^l^") - [ div(«^i), 

Jv 

which cancels out with the same term VmJj^^ in the momentum equation. □ 

5.3.2 The discontinuous Galerkin method 

If the advection term u ■ "Vcr is discretized by the discontinuous Galerkin method, the system 
writes: 



fc=i "-f^fc 



- div v + g div ul+^ + (1 - e) V<+i : Vv + ^^tth (e^^'^') : Vv 

, ^ (^„+r^^ (^„+r^ _ , ^ ^^3^ 



^ / K.n|K(tAri)l:</,^ 
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Proposition l5.2l still holds for solutfons of the system ([75]) . The proof foUows that of the previous 
Section [5.3.11 except for the treament of the jump term, which follows that of the Section l4.1.2l fsee 
also ITT^ . because tth (^A^^^) e (Po)^^^ and 7r,i (e^J-'^'j = is also in (Po)^^^. 

5.4 Other finite elements for {uh,Ph) 

In this section, we review the modifications that apply to the systems in the two previous Sec- 
tions [STU and [^31 when the different mixed finite element spaces for {uh,Ph) proposed in Section r4.3l 
are used instead of Scott- Vogelius. Notice that the conclusions of Table[l]about the conditions that 
the velocity field has to satisfy still hold for the two previous Sections 15.21 and 15.31 with piecewise 
linear approximations of cr/j, ■0,^. 

Other finite elements space for {uh,Ph) than Scott- Vogelius and adequate projections of the 
velocity field (see summary in Table ^ have to be combined with interpolations of the stress field 
(Th,iph using TTh (see the two previous Sections [5.21 and 1^31 above). We give a summary of the tt/i 
projections that are then required in Table [H 

5.4.1 Alternative mixed finite element space for {ufi,Ph) with inf-sup condition 

The situation is very similar to that in Section 14.3.21 Among the mixed finite element space that 
satisfy the inf-sup condition, let us first choose the Taylor-Hood (P2)'' x Pi. Again, because the 
velocity is not even weakly incompressible in the sense of the equation (|55|) . we need to use the 
projection of the velocity field onto the solenoidal vector fields for the treament of some terms in 
the variational formulations. When the advection terms u ■ ^ cr and u ■ V-i/^ are discretized using 
the characteristic method, we define the flow with P™* (m)J) like in ([55|l and use the same systems 
([65]) and ((70|) than above. When the advection terms u ■ 'Vcr and u ■ Wxjj are discretized using the 
discontinuous Galerkin method, we use systems similar to (j65|) and (j70p above, where the jump 
term rewrites (in the conformation-tensor formulation): 

+J2 / iprK)-^iKK+^)i:0+. 

And one still needs to add the so-called Temam correction term ([59|) to the weak formulation. 

We can also use the Crouzeix-Raviart finite elements for velocity (see (j56l) ): {uh,Ph,crh) in 
^CRY X Pq X (Pi, disc) * 2^ ' . Similarly to the advection terms u ■ 'Vcr and u ■ '^'4>, the advection 
term u ■ Vit in the Navier-Stokes equations should then be discretized either using a characteristic 
method with the fiow defined in ([HO]) with any of the projections Ph introduced above for the 
velocity field, or using the discontinuous Galerkin method formulated in the equation ([62|) . 

It is noticeable that choosing the mixed finite elements of Crouzeix-Raviart simplifies all the 
variational formulations presented above in the present Section [51 Indeed, since Vm G (Pq)'^^'' 
and we have the Lemma l5.H it is then unnecessary to project the velocity except in the advection 
terms. For instance, for the conformation-tensor formulation using the discontinuous Galerkin 
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method, the formulation writes: 



-pI+^ divv + qdivu]l+^ + (1 - £)V<+1 : Vv + ^crj^'+i : Vv 



At 



Wi^ ^ ^ 



- ((V<+i)^ri + <+i(V<+i)^) : 4> (74) 



+ ^ / IP, K) . n| [tt, : 0+ + Re |P, «) • n\ K+'l ' {^} • 

Note that the second term in the sum of integrals over edges Ej is due to the use of the Crouzeix- 
Raviart element, and is uncorrelated to the treatment of the advection by a discontinuous Galerkin 
method. 

The discrete free energy estimate (|57| holds. Its proof combines arguments of the proofs above, 
except for the treatment of the upper-convective term in (|74p . This term writes, on any element 
Kk of the mesh (using Lemma \5A\ the fact that Vtt e (Pq)''^'^ and the incomprcssibility ([M]) ): 

/ T7„ /T „ f—n+l\^^\ / _Ji+l T-r„ / _ r7„ "+1 

/ ^^h (^h ■ {I - (cTft ) ) = / CTft : - '^h ■ (cr^ ) Vu^ , 

JKk JKk J-D 

= / : V<+i - / div< 



,n+l 
''h ' 



which vanishes after combination with the extra-stress term in the momentum equation, the latter 
satisfying: 



_ri+l T-7„ n+1 / _ /'_n+l^ tt„ n+1 

Kk JKk 
dxd 



because of the fact that Vm e (Pq) and using Lemma [57T] 

5.4.2 Alternative mixed finite element space for {uh,Ph) witiiout inf-sup 

It is also possible to use finite element spaces for {uh,Ph) that do not satisfy the inf-sup condition 
like in Section 14.3.31 while the stress field is discretized using discontinuous piecewise linear ap- 
proximations. The construction of systems of equations and the derivation of discrete free energy 
estimates then directly follow from the combination of results from the Section 14.3.31 with those 
used above in Section [51 after upgrading the degree of the polynomial approximations for the stress 
field. 
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If we consider the mixed finite element space (Pi)'^ x Pq for {uh,Ph), and if the term u ■ 'Vcr is 
discretized with tlic characteristic method, the system then writes: 



k=l 



-/^'+Miv'i; + <7divw;j+i + (l-e)VM;j+i : '^^ + ■ 



Ne „ 



with a flow X" computed with the projected field P™* (itJJ) through It is noteworthy that, 

for the same reason than above in the equation ()74p . the projection operator tt^i is needed only for 
the discretization of the advection term u ■ ^cr. 

If we consider the mixed finite element space (Pi)'^ x Pi for {uii,ph), it is straightforward to 
rewrite the system (jM)) where the stress field was only piecewise constant, while using the same 
argument as above to see that only the advection term for the stress field needs a projected velocity. 

Remark 5.3. We were not able to retrieve discrete free energy estimates without interpolating 
some terms in the formulations above thanks to the operator tt/j. This operator projects the stress 
(Til (or tpfi) onto (Pq) * 2^ ' . Thus, for the formulations we have considered in this section, the 
interest of using larger dimensional spaces for cr/j ( or rpfi) than (Pq) 2 is not clear. Our aim in 
this section is simply to exhibit discrete formulations with piecewise linear approximations of the 
stress, for which we are able to derive a free energy estimate. 



6 Positivity, free energy estimate and the long-time issue 

Notice that both Propositions 13.11 and 13.21 impose a limitation on the time step which depends 
on the time iteration: < At < cq, where cq = co(w5J,crJJ) is function of a time-dependent data. 
Thus, these existence results are weak insofar as the long-time existence of the discrete solutions 
is not insured, i.e. if J2'^=o "^ol^^, cr^) < 00. 

Yet, for the discretizations introduced above, we have also shown that at each time step, the 
solutions of those discretizations satisfy free energy estimates. This will now allow us to prove the 
long-time existence of the discrete solutions: 

Proposition 6.1. For any initial condition (it^,cr^) with cr^ symmetric positive definite, there 
exists a constant ci = ci (M^,cr^) > such that, for any time step < At < ci, there exists, for 
all iterations n g N, {u^'^^ , cr^'^^) which is the unique solution to the system p7p (resp. (I39p ) with 
(T^^^ symmetric positive definite. 
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Vu. 



in (Po)'^' 



not in (Po)''' 



TT/i (crJJ) in time derivative 
(incl. flux term in DG) 



l.disc 



TTh i'tph.) in time derivative 
(incl. flux term in DG) 
+ implicit source term nh (e^'^h ' ) in OB 



NS 



implicit coupling term tt^ ^e'^'*^ ^ 



TTfi (crJJ) in time derivative 
(incl. flux term in DG) 

+ implicit coupling terms 
( TTH {(t]1+') in NS, OB ) 



TTh (tph) in time derivative 
(incl. flux term in DG) 
+ implicit source term tt^, (^e~^i' 
+ implicit coupling terms 
( (e'^^^') in NS, tt^ {xPl+^) in OB ) 



OB 



Table 3: Summary of projected terms in the Navier-Stokes (NS) and Oldroyd-B (OB) equations 



for CTh/iph in (I 



l^disc 



Proof of Proposition 16. Jl Like in the proof of Proposition 13.11 we will proceed with the proof for 
system ([57| only, using its rewriting as system pSj) . 

The proof is by induction on the time index n. With the notation of the proof of Proposition l3.H 
for a fixed n = . . . iV^ — 1 and for a fixed vector Y" of values in the subset S*^ of jjZTVD+sAf/f 
(standing for the nodal and elementwise values of a field (mJJ,(tJJ) with crJJ symmetric positive 
definite), we define like in the proof of Proposition 13.11 (using the implicit function theorem) a 
function R : At G [0,co) i?(At) e rstVo+stVk (where cq = co«,(t'^')) such that: 

VAi e [0,co),Q(At,i?(Ai)) = 0, 

where Q is deflned by (gB]). For any At e [0,co), R{At) e stands for the nodal and 

elementwise values of a field {uh{At),crh(At)) (with (T/j(At) symmetric positive definite) that is 
solution to the system P5|l . 

Then, by Proposition 14. II the solution (u/i(A<), cr/j(A<)) satisfies a free energy estimate: 

F(«,(Ai), <T,(Ai)) < F«, ^',') . (76) 

Using the fact that all norms are equivalent in the finite-dimensional vector space ^"^^o+sNk ^ 
and that, for Q < < 1 — i, we have i/ x < x — ln(x) ,Vx > 0, we obtain that there exists two 
constants a > and (3 > (independent of At), such that: 

a||i?(At)|| < FiuhiAt), cThiAt)) + (3 . (77) 

Let us define the function D: 

D:Ate [0, Co) — > B{Y") + A{R{At)) + (VyA)i?(At) e R2^«+3A'ic_ 
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We recall that (see (gT])), with the new notations: VyQi^t, R{At)) = I + AtD{At). Us- 
ing ([75]). ([77)1 and the fact that the discrete free energy is non- increasing, the function D satisfies: 

\\DiAt)\\<\\B\\\\Y^ + {\\A\\ + \\VYA\\)\\RiAt)\\, 

<{\\B\\ + \\A\\ + \\WyA\\)-(F{uI,ctI)+P), 

a 

< (IIBII + P!| + II Vy^||)i (fK, ^"J -f /3) . 

a 

This shows that there exists a constant ci = ci , cr°) > such that, for any time step < At < 
ci, the matrix VyQiAt, R{At)) is invertible. Using the implicit function theorem, this implies 
that, for any time step < At < ci, the system admits a solution (mJJ^^, cr,"^^) with crj^^^ 
symmetric positive definite at all iterations n £ N. □ 

A similar result can be proven for the log- formulations (|4ip and (|^^ : 

Proposition 6.2. For any initial condition (m^j,i/'^J, there exists a constant ci = ci (m^j,!/?^) > 
such that, for any time step < At < ci, there exists, for all iterations n G N, (m^^^, i/'J^''"^) which 
is the unique solution to the system (I4ip (resp. (|42p ). 

Proof of Proposition 1 6'. The proof of Proposition 16.21 is fully similar to that of Proposition 16.11 
using for Q{At, Y) and D{At) slightly modified expressions as explained for the proof of Propo- 
sition [3T2l The entropic term in the free energy still helps at bounding the norm of the vector of 
nodal-elementwise values for {uh,,ipii) like in (|77p using the following scalar inequality, which is 
true for any fixed e (0, 1]: Vx G M, — x + 1 > i/|x|. □ 

From Propositions 16.11 and 16. 2[ we have the global-in-time existence of solutions to those dis- 
cretizations of the Oldroyd-B system presented above which satisfy a discrete free energy estimate. 

The log-formulation actually also satisfies the following long-time existence result, using the 
fact that the a priori estimates can be obtained without requiring the stress tensor field to be 
positive definite: 

Proposition 6.3. For any initial condition (u^, ^A/l), and for any constant time step At > 0, there 
exists, for all iterations n CzN, {u^^^^ which is solution to the system ()4ip (resp. ()42p ). 

Proof of Proposition \6.3[ We will proceed with the proof for system (|4ip only, using its rewriting 
as system (|45|) . Note already that, since the derivation of discrete free energy estimates for the 
system (|^IT) does not require the solution tp]^^^ and the test function to be non-negative like in 
the derivation of discrete free energy estimates for the system ((S7)) . then the manipulations used to 
derive the free energy estimate ()53p can also be done a priori for any function in the finite element 
space. 

Let us consider a fixed time index n and a given couple {uf^,xp'^) e {¥2)^^=0 ^ i^o) We 
equip the Hilbert space 

(IP2)div=o ^ (IPo)~^ with the following inner product: 

((^'l>0l); (^2,02)) = / Vi ■ V2 + 4>i : , 

Jv 
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for all {vi,4>i),{v2,4>2) £ (IP'2)div=o ^ (^o) ^ , and denote by || • || the associated norm. Let 
us introduce the mapping T : (P2)div=o x (Po)^^^ ^ {^2)^,^=0 x (Po)^^^ defined by: for 
iu,'iP) e (P2)L=o X in)^, for any test function e (P2)^iv=o x i^a)^ , 

{T{u, i>); {v, 0)) - y" Re + < • Vm^ • t; + (1 - e) Vm : Vt; + :^e^ : Vv 



At 

where 17 and B are associated with the decomposition of the velocity gradient Vtt as explained in 
Lemma 12.21 

Vu = n + B + Ne-'f'. 
If {ul+\xp"+^) is solution to (051), then we have: for all (t>,</.) e (P2)div=o ^ (Po)^^^, 

{H<\^l^'yAv,ct>))=0. (78) 

Considering and using 2Wi (^''^ ~ ■^)) ^^st functions, we get the following inequality 

after similar manipulations to those in the proof of Proposition [ 



Re /■ , ,9 e f , , Re 



> ^ \v 



2Wi ., 2Wi 



' tr(e<^ -^)-^ \K? - — / tr(e^^ - ^n) (79) 



Jv ^ Jt> 



V 

r'l' + ^tr (e<^ + e-<^ - 2/) 



Let us now assume that, for any a > 0, the mapping T has no zero (u^^'^ ^ V'/J^^) satisfying ((75)) 
in the ball 

6„ - {(t;,0) e (P2)iv=o X (Po)"^, 11(^^,0)11 < a} ■ 

Then, for such a real number a > to be fixed later, we define the following continuous mapping 
from Ba onto itself: 

Q{v,c^) = -a^g^,V(^,0) e (P2)^i..o X (Po)^. 
ll-^(«,0)ll 

By the Brouwer fixed point theorem, Q has a fixed point in Ba- Let us still denote that fixed point 
(v, 0) for the sake of simplicity. By definition, it satisfies: 

g(y, 0) = (t;, 4>) e 6a and \\g{v, 0)|| a. (80) 
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Using the equation ([50]) with the inequahty ([7^ . we have: 



a V' V ' 2Wi 



a 



+ ^ - <p + At ^(1 - £)|V<+ip + ^ tr (e<^ + e-'f - 2l) . (81) 

Using the scalar inequality — x > |a;|, Va; g R, we have: 

/ tr(e<^-0 + /)>^ / |A,|,V0e (Po)™, (82) 

J-D ^ Jv 



where (Ai)i<i<d are functions depending on <p such that, for all x ^ V, (Ai(a;))i<j<d are the d 
(non-necessarily distinct) real eigenvalues of the symmetric matrix <f>{x). 

Now, since (IP2)div=o ^ (Po) ' ^ ' is finite-dimensional, all norms are equivalent. So there exist 
71,72 > such that, for all (-*;,</.) e (P2):li,=o x (Po)"^: 

V 



7111(^,0)11 <(/ l^^r) +11 max |A,(a;)|||,,<72||(t;,0)||, (83) 



where it is easy to prove that || maxi<i<rf |Ai(a;)| ||oo defines a norm in the vector space L°° (V, 5(M''^'^)) . 
Using the equation ((5^ with the norm equivalence we obtain: 



Re 

T 



. ,'Re e 1 
> mm — , ■ 



2 '2Wi||maxi<,<d|A,(a;)||| 



(84) 

^I^P + II max |A.(a:)||UE/,l^^l) ' (85) 



>min('5£,-£-^_ \^^^^,„ ] (/^H^ + E /^|A.p) , (86) 



2 ' 2Wi II maxi<j<<i |Ai(a;)||| 

Re ^ 1 

T'2Wi||maxi<,<d|A,(a;)||| 



ll(^,0)f, (87) 



-"""'T'2Wi^J" 



If we now choose a large enough so that: 
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we get: 

^ [ \v\^ + — [ tr(p> - + J) _ / \ul\^-^ [ tY{e^i -iPl + I) 
+ J^^\^- K\' + At y^(l - e)| V<+ip + ^ tr (e-^ + e'^ - 2l) > 0. 

This is obviously in contradiction with since, for aU 4> G (Pq) * ^'^ ' , we have tr((^e'^ — 0) > 
in virtue of the scalar inequality x{e^ — 1) > 0, Vx G M. 

Thus, for any At > 0, if we choose a sufhciently large, the mapping T has a zero {u^^^^ 
satisfying ((781) in the ball Ba- □ 



Notice that Proposition 16.31 does not ensure the uniqueness of solutions. The fact that we are 
able to prove such a stability result without any assumption on the timestep for the log-formulation, 
and not for the classical formulation, may be related to the fact that the log-formulation has been 
reported to be more stable than the classical formulation (see [23]). 

Remark 6.1 (Other positivity preserving schemes). There exist other means than using the log- 
formulation to preserve the non-negativity of the conformation tensor. A very natural way of 
preserving the non-negativity is to use a Lie-formulation of the upper convective derivative term. 
For instance, following ]32fj . one can write the following discretization of ([32| using Scott- Vogelius 
elements for (uh,Ph) and piecewise constant approximations for (Th'- 

0= / Re['^^!^^^^+ul-Vul+^\-v 



V 



( -{I- Atn, (V<+^))-^ « o X-jt-)) [I Atn, (V<'+^))-^ 

At 



(90) 





Wi^ ' 

where the function X"(t) is defined by psp . The system ()90p admits a solution such that (J — 
Atnh (^■'^/J^^))^^ is well-defined, provided At is sufficiently small (hut possible very small when 
II VmJj'^^II is large). Besides, taking cf) as the characteristic function of some element Kk, we have 
the following equality inside Kk : 

^ + Wi) "^''^^ = " c^^^'))"' ^< ° (-^ - c^^r'))""^ + (91) 

Then it is clear that the system (|90p preserves the non-negativity of cr^f^. Moreover, it is possible to 
derive the free energy estimate (|49p for the system (j90p . It suffices to take as a test function for 
the stress: 
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and to proceed to the derivation of a free energy estimate using both ideas of the present work and 
of the work 1321, after noting that: 

tr (vr. (V«r^)" (I ~ (-r^)-^) (V<+^) ((l + <Tr^ - 1^/)) > 0, (92) 

the proof of which is fully similar to the proof of @, using the fact that (l + cr^^^ — -^J is sym- 
metric positive definite (provided H^t is sufficiently small) andir^ (VitJJ^"'^)"^ — i'^'h^^)^^^ '^h (^'^Ji 
is symmetric positive semi- definite. 

A Some properties of symmetric positive definite matrices 
A.l Proof of Lemma 11.11 

Formula ([6]), ([7]) and ([8]) are simply obtained by diagonalizing the symmetric positive definite 
matrix cr, and using the inequalities: Va;, y > 0, ln(a;?;) = Inx + In?/, a; — 1 > Ina; and x -\-\/x>2. 

Let us now prove Formula By diagonlization, we have cr — Q,^ Dil with fl orthogonal and 
D diagonal strictly positive, which gives: 

tr(crr"^) = tr(Sl^\/DVI?OT"i), (93) 
= triVDriT-^n^VD), (94) 
> 0, (95) 

because ^/DflT~^^'^ ^/D is clearly a symmetric positive definite matrix. 
Likewise, we have: 

det(cn--i) = detin^ DflT-^), 

And noting that A — ^/DflT^^ri^^/D is symmetric positive definite, the proof of (fTU|) is then 
equivalent to show: 

ln(det(yl)) <tr(A- J), 

for any symmetric positive definite A, which derives from ^ and ([7]). 

It remains to prove Formula By diagonalization, we write cr = DO and r = R^AR 

with and R orthogonal, and D and A diagonal strictly positive. Let us introduce the orthogonal 
matrix fl — OR^ . We denote by Di (resp. A^) is the («, i)-th entry of D (resp. of A). 
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Let us first prove Formula PT|) in dimension 2. We set (p = fif^ {(f> £ [0, 1]), so that — i^~'P)^ 
^21 = (1 ^ 'A) ^'^d — 4>- We have 

J((/)) = tr ((In cr — In r) cr — (T + r) , 

= ^ A In A - A + A, - ^ nf^D, In Aj, 

= ^ A hi A - A + A, - (/.^ A In A, - (1 - ^ A hi A^-. 

The functional J is an afhne function defined on the interval [0, 1]. For = 0, we have 

J(0) = ^ D, In A - Di + - ^ A hi Aj , 

^i:A(%f-i-.n(4))^o, 

for s the permutation of indices 1^2. We have used the scalar inequality a; — 1 > Inx for all 
X > 0. Likewise, for 0=1, 

J(l) = ^ D, In -D,+K,-Y, In A„ 

The functional J is thus non negative on [0, 1], which ends the proof in dimension 2. 
More generally, in dimension d, let us define the following affinc functional: 

J({f7^^.})=tr((ln^-lnr)<T-(^-T)), 

= ^ A In A - A + A,; - A In A,. 

The functional J is defined on the convex polyhedron V = G M.f , Vi J2j ^Ij = 1> Vj Y.^ ^% = ^ 

and Q = e R'^ , 5757-'" = j| is a subset of T'. All maxima and minima of J on P are reached 

at vertices of the convex polyhedron V (Theorem 10.3-2 in lOj). Besides, {57?^ } is a vertex (that 
is, an extremal point) of P if, and only if, for all 1 < i < d, there exists one, and only one, j = s(i), 
1 < j < c^, such that f^^^fj^ 7^ (Theorem 10.3-1 in [IH])- Thus, the vertices of V are exactly the 
{17?.} such 



02 _ J 1' J = S(«) 



, 0, s(«) 

where s is a permutation on {1, . . . , d}. Now, for any permutation s, we have: 

J = J2^i In D, - D, + A, - Di In A^(,) 

= EA(^-i-..(%f^))^o, 
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using the scalar inequality a; — 1 < In a; for all a; > 0. Hence the result: J({ilfj }) > 0, V{r2fj } e Q. 
A. 2 Proof of Lemma Fl.21 

Since cr e C^([0,T)) is symmetric positive definite, it can be decomposed as: Vi e [0,T) 

cr(t) = R{tfK{t)R{t) 

where the unitary matrix R and the diagonal matrix A are continuously differentiable. Let us 
compute the following derivatives: 

d , dR^ , , ^ dA . 1 -r, ,dR 

— Incr = — — lnAi? + i?^^— A"ii? + i?^lnA— , 
at at dt dt 

d dR^ . dA T , di? 

— cr = ——AR + R^—R + R^A—. 
dt dt dt dt 



Thus, we obtain: 



and 



Now, for all times t, 



R^R = I, 

hence (IT^ . The proof of is similar. 

B Proof of Lemma 12.11 

Let us introduce 

to — inf{t > 0, cr is not symmetric positive definite}, 

with convention to = oo if {t > 0, cr is not symmetric positive definite} = 0. Since cr(t = 0) is 
symmetric positive definite, it remains so at least for small times < t < At, by continuity of the 
eigenvalues with respect to the time variable t. Thus, to > At > 0. 

Let us assume that to < oo. First, one can define the logarithm Incr of the matrix cr, which 
satisfies the equation for in system ([23| for t G [0, to). Taking the trace of the equation for ip in 
system (^5]) . we get for Incr: 

— In deter = —tr(cr"i -/), (96) 
I^t Wi ^ ^ ^ ' 

where we have introduced the convective derivative 

D-t=[dt+^^-^^ 
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Besides, for any positive definite matrix cr ^, we have: 

^^^^^> (det<T-i)i/^ (97) 
d 

which foUows from the convex inequahty between geometrical and arithmetical means. Thus, 
combining and ((571) . we get, on the time interval [0,to): 

^(dct<T)i/'^ = i(det<T)i/'^^lndet^ ^ Wi " (det<T)V'') . (98) 

Now, by continuity of the eigenvalues with respect to t, one eigenvalue at least converges to zero 
as t Iq , which implies det cr — s- 0"*". Then, there exists 77 > such that, for times to ^ rj < t < to, 
we have: 

< deter < 1, 

and by ^ 

-^(det cr)!/'^ > 0. (99) 

But then, tg cannot be the first time when det cr = 0, otherwise one should have -^(det (Ty/'^{t^) < 
0, which contradicts Thus to = 00 which ends the proof of Lemma \TT[ 



C Proof of Lemmas [23] and [33] 

Let <T be any piecewise constant function with symmetric positive definite matrix values cr G 
S'^{{PoY^'^). We define the following linear applications for functions N with skew-symmetric real 
matrix values TV e A{R'^'^''): 

f^:N^^{Na-^-<T-^N) (100) 
g^:N ^^{Na-' +(T-'N). (101) 

The linear application go- is an automorphism in A(M.'^^'^) and the linear application /„■ is a 
one-to-one onto function from ^(M'^^'*) to the set of symmetric real matrices with null trace 
cr e 5"((Po)'^'''').Both applications are of rank 1. 

Now, for any x Cz V, the symmetric positive definite matrix cr(a;) can be diagonalized with an 
orthogonal matrix 0{x) E 0{M.'^^'^) {0^0 = I = 00^) and a diagonal positive definite matrix 
D{x) e D((M;)''): 

ct{x) = 0{xfD(x)0{x). 

So, if we change the basis, we can equivalently work with a diagonal positive definite matrix 
cr(a;) e D((R;)'^). 

Let u'f^'^^ be some vector field in (Pfe)^ with fc e N and ijj^^^ G S ^(Po) * ' ) elementwise 
constant symmetric rank-two-tensor field. The rank-two-tensor field Vtt^^^ takes its values in 
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(Pfc_i^disc)'^^'*- For X £ V, Vm^'^^(cc) decomposes elementwise in the basis where the matrix 



-h 

d(d+l) 



of the hnear apphcation associated with the symmetric matrix V'/J^^ G ^(Po) ^ is diagonal. 
More precisely, in this basis, ^u^J^^{x) uniquely decomposes as a direct sum in the supplementary 
spaces (of respective dimension 1, 2 and 1): 

This decomposition reads: 

+ di.g(V<«(.)) + ^""''-'-7""'W^ (103) 

Since e~'^h* is an elementwise constant diagonal positive definite matrix, in order to get the 
following decomposition at each x ^V: 

Vul+\x) ^ ni+\x)+Bl+\x)+Nl+\x)e-^l^\x) (104) 

BT\x) (105) 



+ ni+\x) + tlJ^^B}^ ±^ (log) 

we need to do the following identification: 

Bl+\x) = di^g{Vul+\x)), (107) 

Nl-^'i^)-!-},.,.,^, + _diag(Vnr^(.))) , (108) 



(-) \^ 2 

which is unique because of the requirement that B^^'^^{x) commutes with the diagonal matrix 

With the choice above, the matrices fij^^^, BJ^^^, ATJ^"*"^ have entries in (Pfc-i.disc)^^^ because 
g-V-fc+^la:) (Po)^t'^. Moreover, we clearly have: tr(B5J;+"^) = divM5J+\ Besides, notice that the 
triplet (JIJ^'*"^, BJJ'*'^, TVJJ^^) is unique provided the change of basis for e^^i^ to be diagonal is 
unique, which means provided that e^'^^^ is not singular (all the eigenvalues should be distinct 
one another). 



D Proof of Lemma 15.1 



d(d+i) 

For any G iodise) ^ , we want to show that: 



(t>h -(Ph^ / i4>h) ■■ 4)1,, V^^ e (Po) ^ . 
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It is enough to prove Lemma |5. II on each simplex Kk G Tii and in the scalar case. Let {xi)i<i<3 
be the vertices of the simplex Kk and {'4'i)i<i<3 the corresponding (linear) basis functions in Pi. 
Then, the function ^/J/^-j. G Pi reads 

0ft|xfc(a;) = cf>f^{x i)tpi{x) + cf>f^{x2)tp2{x) + 4)h{x3)tp3{x), Vx G Kk- 

For every (f)f^ £ Pq, 

because J^^ ipi = Moreover, 0/i|/^fc S Pi, hence 

which means 

JKk Jh'k JKk 

Acknowledgements We would like to thank A. Ern and J.W. Barrett for fruitful discussions. 

References 

[1] D.N. Arnold and J. Qin. Quadratic velocity/linear pressure Stokes elements. In R. Vichn- 
evetsky and R.S. Steplemen, editors, Advances in Computer Methods for Partial Differential 
Equations, volume VII, 1992. 

[2] J. Baranger and A. Machmoum. Existence of approximate solutions and error bounds for 
viscoelastic fluid flow : Characteristics method. Comput. Methods Appl. Mech. Engrg., 148:39- 
52, 1997. 

[3] J.W. Barrett and S. Boyaval. Global solutions to some regularized free-energy-dissipative 
Oldroyd-like models. In progress. 

[4] J.W. Barrett, C. Schwab, and E. Suli. Existence of global weak solutions for some polymeric 
flow models. Math. Models and Methods in Applied Sciences, 15(6):939-983, 2005. 

[5] A.N. Beris and B.J. Edwards. Thermodynamics of flowing systems with internal microstruc- 
ture. Oxford University Press, 1994. 

[6] J. Bonvin, M. Picasso, and R. Stenberg. GLS and EVSS methods for a three fields Stokes 
problem arising from viscoelastic flows. Comp. Meth. Appl. Mech. Eng., 190:3893-3914,2001. 

[7] F. Brezzi, J.R. Douglas, and L.D. Marini. Two families of mixed finite elements for second 
order elliptic problems. Numer. Math., 47:217-235, 1985. 



RR n° 6413 



50 



Doijaral cl at 



[8] F. Brezzi, J.R. Douglas, and L.D. Marini. Recent results on mixed finite element methods for 
second order elliptic problems. In A.V. Balakrishnan, A. A. Dorodnitsyn, and J.L. Lions, edi- 
tors. Vistas in Applied Mathematics: Numerical Analysis, Atmospheric Sciences, Immunology, 
pages 25-43, 1986. 

[9] F. Brezzi and J. Pitkaranta. On the stabilization of finite element approximations of the 
Stokes equations. In W. Hackbusch, editor. Efficient Solution of Elliptic System, pages 11-19, 
1984. 

[10] P.G. Ciarlet. Introduction a I'analyse numerique matricielle et d V optimisation. Masson, 1985. 
In French. 

[11] R. Codina. Comparison of some finite element methods for solving the diffusion-convection- 
reaction equation. Comp. Meth. Appl. Meek. Engrg., 156:185-210, 1998. 

[12] M. Crouzeix and P.A. Raviart. Conforming and non-conforming finite element methods for 
solving the stationary Stokes equations. RAIRO Anal. Numer., 3:33 75, 1973. 

[13] A. Ern and J.-L. Guermond. Theory and practice of finite elements. Springer Verlag, New- 
York, 2004. 

[14] A. Fattal, O.H. Hald, G. Katriel, and R. Kupferman. Global stability of equilibrimii manifolds, 
and "peaking" behavior in quadratic differential systems related to viscoelastic models. J. Non- 
Newtonian Fluid Mech., 144:30-41, 2007. 

[15] R. Fattal and R. Kupferman. Constitutive laws for the matrix-logarithm of the conformation 
tensor. J. Non-Newtonian Fluid Mech., 123:281-285, 2004. 

[16] R. Fattal and R. Kupferman. Time-dependent simulation of viscoelastic flows at high Weis- 
senberg number using the log-conformation representation. J. Non-Newtonian Fluid Mech., 
126:23-37, 2005. 

[17] E. Fernandez-Cara, F. Guillen, and R.R. Ortega. Mathematical modeling and analysis of 
viscoelastic fluids of the Oldroyd kind. In P. G. et al. Ciarlet, editor. Handbook of numerical 
analysis. Vol. 8: Solution of equations in M" (Part 4-)- Techniques of scientific computing 
(Part 4)- Numerical methods of fluids (Part 2)., pages 543-661. Elsevier, 2002. 

[18] J.L. Guermond. Stabilization of galcrkin approximations of transport equations by subgrid 
modeling. Math. Model. Numer. Anal., 33(6):1293 1316, 1999. 

[19] C. Guillope and J.C. Saut. Existence results for the flow of viscoelastic fluids with a differential 
constitutive law. Nonlinear Analysis, Theory, Methods & Appl., 15(9):849-869, 1990. 

[20] D. Hu and T. Lelievrc. New entropy estimates for the Oldroyd-B model, and related models, 
2007. Available at http://arxiv.org/abs/inath.NA/0703198, to appear in Commun. Math. 
Sci. 



INRIA 



Dissipative Oldroyd-B schemes 



51 



[21] T.J.R. Hughes and L.P. Franca. A new finite element formulation for CFD: VII the Stokes 
problem with various well-posed boundary conditions: Symmetric formulations that converge 
for all velocity/pressure spaces. Comp. Meth. App. Mech. Eng., 65:85-96, 1987. 

[22] M.A. Hulsen. A sufhcient condition for a positive definite configuration tensor in differential 
models. J. Non-Newtonian Fluid Mech., 38:93-100, 1990. 

[23] M.A. Hulsen, R. Fattal, and R. Kupferman. Flow of viscoelastic fluids past a cylinder at 
high Weissenberg number: stabilized simulations using matrix logarithms. Journal of Non- 
Newtonian Fluid Mechanics, 127(l):27-39, 2005. 

[24] B. Jourdain, C. Le Bris, T. Lelievre, and F. Otto. Long-time asymptotics of a multiscale 
model for polymeric fluid flows. Archive for Rational Mechanics and Analysis, 181(1):97-148, 
2006. 

[25] N. Kechkar and D. Silvester. Analysis of locally stabilized mixed finite element methods for 
the Stokes problem. Mathematics of Computation, 58(197):1-10, 1992. 

[26] R.A. Keiller. Numerical instability of time-dependent flows. J. Non-Newtonian Fluid Mech., 
43:229-246, 1992. 

[27] R. Keunings. Simulation of viscoelastic fluid flow. In C. Tucker, editor, Fundamentals of 
Computer Modeling for Polymer Processing, pages 402-470. Hanse, 1989. 

[28] R. Keunings. A survey of computational rheology. In D.M. Binding et al., editor, Proc. 13th 
Int. Congr. on Rheology, pages 7-14. British Society of Rheology, 2000. 

[29] R. Kupferman, C. Mangoubi, and E. Titi. A Beale-Kato-Majda breakdown criterion for an 
Oldroyd-B fluid in the creeping flow regime. To appear in Comm. Math. Sci. Preprint available 
at ht tp : //arxiv . org/PS_cache/arxiv/pdf /0709/0709 . 1455vl . pdf . 

[30] Y. Kwon. Finite element analysis of planar 4:1 contraction flow with the tensor-logarithmic for- 
mulation of differential constitutive equations. Korea- Australia Rheology Journal, 16(4):183- 
191, 2004. 

[31] Y. Kwon and A.V. Leonov. Stability constraints in the formulation of viscoelastic constitutive 
equations. J. Non-Newtonian Fluid Mech., 58(l):25-46, 1995. 

[32] Y. Lee and J. Xu. New formulations positivity preserving discretizations and stability analysis 
for non-Newtonian flow models. Comput. Methods Appl. Mech. Engrg., 195:1180-1206, 2006. 

[33] A.I. Leonov. Analysis of simple constitutive equations for viscoelastic liquids. J. non-newton. 
fluid mech., 42(3):323-350, 1992. 

[34] F.-H. Lin, C. Liu, and P.W. Zhang. On hydrodynamics of viscoelastic fluids. Comm. Pure 
Appl. Math., 58(11):1437-1471, 2005. 



RR n° 6413 



52 



Doijaral cl at 



[35] P.L. Lions and N. Masmoudi. Global solutions for some Oldroyd models of non-Newtonian 
flows. Chin. Ann. Math., Ser. B, 21(2):131-146, 2000. 

[36] A. Lozinski and R.G. Owens. An energy estimate for the Oldroyd-B model: theory and 
applications. J. Non-Newtonian Fluid Meek., 112:161-176, 2003. 



[38] O. Pironneau. On the transport-difFusion algorithm and its application to the Navier-Stokes 
equations. Numer. Math., 3:309-332, 1982. 

[39] D. Sandri. Non integrable extra stress tensor solution for a flow in a bounded domain of an 
Oldroyd fluid. Acta Mech., 135(l-2):95-99, 1999. 

[40] L.R. Scott and M. Vogelius. Norm estimates for a maximal right inverse of the divergence 
operator in spaces of piecewise polynomials. RAIRO Model. Math. Anal. Numer., 19:111-143, 



[41] R. Temam. Sur I'approximation des equations de Navier-Stokes. C. R. Acad. Sc. Paris, Serie 
A, 262:219-221, 1966. 

[42] P. Wapperom and M.A. Hulsen. Thermodynamics of viscoelastic fluids: the temperature 
equation. J. Rheol, 42(5):999-1019, 1998. 

[43] P. Wapperom, R. Kcunings, and V. Lcgat. The Backward-tracking Lagrangian Particle 
Method for transient viscoelastic flows. J. Non-Newtonian Fluid Mech., 91:273-295, 2000. 




Beyond Equilibrium Thermodynamics. Wiley, 2005. 



1985. 



INRIA 




Unite de recherche INRIA Rocquencourt 
Domaine de Voluceau - Rocquencourt - BP 105 - 78153 Le Chesnay Cedex (France) 

Unite de recherche INRIA Futurs : Pare Club Orsay Universite - ZAC des Vignes 
4, rue Jacques Monod - 91893 ORSAY Cedex (France) 
Unite de recherche INRIA Lorraine : LORIA, Technopole de Nancy-Brabois - Campus scientifique 
615, rue du Jardin Botanique - BP 101 - 54602 Villers-les-Nancy Cedex (France) 
Unite de recherche INRIA Rennes : IRISA, Campus universitaire de Beaulieu - 35042 Rennes Cedex (France) 
Unite de recherche INRIA Rhone-Alpes : 655, avenue de I'Europe - 38334 Montbonnot Saint-Ismier (France) 
Unite de recherche INRIA Sophia Antipolis : 2004, route des Lucioles - BP 93 - 06902 Sophia Antipolis Cedex (France) 



Editeur 

INRIA - Domaine de Voluceau - Rocquencourt, BP 105 - 78153 Le Chesnay Cedex (France) 

http://www.inria.fr 
ISSN 0249-6399 



